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Abstract 



Background The primary visual cortex of many mammals contains a continuous representation of visual space, 
with a roughly repetitive aperiodic map of orientation preferences superimposed. It was recently found that 
orientation preference maps (OPMs) obey statistical laws which are apparently invariant among species widely 
separated in eutherian evolution. Here, we examine whether one of the most prominent models for the 
optimization of cortical maps, the elastic net (EN) model, can reproduce this common design. The EN model 
generates representations which optimally trade of stimulus space coverage and map continuity. While this 
model has been used in numerous studies, no analytical results about the precise layout of the predicted OPMs 
have been obtained so far. 

Results We present a mathematical approach to analytically calculate the cortical representations predicted by 
the EN model for the joint mapping of stimulus position and orientation. We find that in all previously studied 
regimes, predicted 0PM layouts are perfectly periodic. An unbiased search through the EN parameter space 
identifies a novel regime of aperiodic OPMs with pinwheel densities lower than found in experiments. In an 
extreme limit, aperiodic OPMs quantitatively resembling experimental observations emerge. Stabilization of 
these layouts results from strong nonlocal interactions rather than from a coverage-continuity-compromise. 
Conclusions Our results demonstrate that optimization models for stimulus representations dominated by 
nonlocal suppressive interactions are in principle capable of correctly predicting the common 0PM design. They 
question that visual cortical feature representations can be explained by a coverage-continuity-compromise. 



Introduction 

The pattern of orientation columns in the primary visual cortex (VI) of carnivores, primates and their 
close relatives are among the most intensely-studied structures in the cerebral cortex and a large body of 
experimental (e.g. p^|l3]) and theoretical work (e.g. p^}|39] ) has been dedicated to uncovering its 
organization principles and the circuit level mechanisms that underlie its development and operation. 
Orientation preference maps (0PM) exhibit a roughly repetitive arrangement of preferred orientations in 
which adjacent columns preferring the same orientation are separated by a typical distance in the 



millimeter range [2j-[5 10 . This range seems to be set by cortical mechanisms both intrinsic to a particular 
area but potentially also involving interactions between different cortical regions |41j. The pattern of 
orientation columns is however not strictly periodic because the precise local arrangement of preferred 
orientation never exactly repeats. Instead, orientation preference maps appear as organized by a spatially 
complex aperiodic array of pinwheel centers, around which columns activated by different stimulus 



orientations are radially arranged like the spokes of a wheel [2||5,10 . The arrangement of these pinwheel 
centers, although spatially irregular, is statistically distinct from a pattern of randomly positioned 
points [38] as well as from patterns of phase singularities in a random pattern of preferred 



orientations 32 , 36 , 38 42 with spatial correlations identical to experimental observations 38 , 42| . This 
suggests that the layout of orientation columns and pinwheels although spatially aperiodic follows a 
definite system of layout rules. Cortical columns can in principle exhibit almost perfectly repetitive order 
as exemplified by ocular dominance (OD) bands in the macaque monkey primary visual cortex |43j[44|. It 
is thus a fundamental question for understanding visual cortical architecture, whether there are layout 
principles that prohibit a spatially exactly periodic organization of orientation columns and instead enforce 
complex arrangements of these columns. 

Recent comparative data has raised the urgency of answering this question and of dissecting what is 
constitutive of such complex layout principles. Kaschube et al. quantitatively compared pinwheel 
arrangements in a large data set from three species widely separated in the evolution of eutherian 



mammals 38 . These authors found that the spatial statistics of pinwheels are surprisingly invariant. In 
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particular, the overall pinwheel density and the variability of pinwheel densities in regions from the scale of 
a single hypercolumn to substantial fractions of the entire primary visual cortex were found to be virtually 
identical. Characterizing pinwheel layout on the scale of individual hypercolumns, they found the 
distributions of nearest-neighbor pinwheel distances to be almost indistinguishable. Further supporting 
common layout rules for orientation columns in carnivores and primates, the spatial configuration of the 
superficial patch system 45 and the responses to drifting grating stimuli were recently found to very 
similar in cat and macaque monkey primary visual cortex 46^. 

From an evolutionary perspective, the occurrence of quantitatively similar layouts for OPMs in primate 
tree shrews and carnivorous species appears highly informative. The evolutionary lineages of these taxa 
diverged more than 65 million years ago during the basal radiation of eutherian mammals [47}|49] . 
According to the fossil record and cladistic reconstructions, their last common ancestors (called the 
boreo-eutherial ancestors) were small-brained, nocturnal, squirrel-like animals of reduced visual abilities 
with a telencephalon containing only a minor neocortical fraction [47||50]. For instance, endocast analysis 
of a representative stem eutherian from the late cretaceous, indicates a total anterior-posterior extent of 
4mm for its entire neocortex [47{[5Q]. Similarly, the tenrec {Echinops telfari)^ one of the closest living 



relatives of the boreoeutherian ancestor 51,52 , has a neocortex of essentially the same size and a visual 
cortex that totals only 2mm^ 47 . Since the neocortex of early mammals was subdivided into several 
cortical areas [47j and orientation hypercolumns measure between 0.4mm^ and 1.4mm^ 38 , it is difficult 
to envision ancestral eutherians with a system of orientation columns. In fact, no extant mammal with a 
visual cortex of such size is known to possess orientation columns [53]. It is therefore conceivable that 
systems of orientation columns independently evolved in laurasiatheria (such as carnivores) and in 
euarchonta (such as tree shrews and primates). Because galagos, tree-shrews, and ferrets strongly differ in 
habitat and ecologically relevant visual behaviors, it is not obvious that the quantitative similarity of 
pinwheel layout rules in their lineages evolved driven by specific functional selection pressures (see 54 for 
an extended discussion). Kaschube and coworkers instead demonstrated that an independent emergence of 
identical layout rules for pinwheels and orientation columns can be explained by mathematically universal 
properties of a wide class of models for neural circuit self-organization. 

According to the self-organization scenario, the common design would result from developmental 
constraints robustly imposed by adopting a particular kind of self-organization mechanism for constructing 
visual cortical circuitry. Even if this scenario is correct, one question still remains: What drove the 
different lineages to adopt a similar self-organization mechanism? As pointed out above, it is not easy to 
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Figure 1: The dimension reduction framework. In the dimension reduction framework, the visual 
cortex is modeled as a two-dimensional sheet that twists in a higher-dimensional stimulus (or feature) space 
to cover it as uniformly as possible while minimizing some measure of continuity (left). In this way, it 
represents a mapping from the cortical surface to the manifold of visual stimulus features such as orientation 
and retinotopy (right). 



conceive that this adoption was favored by the specific demands of their particular visual habitats. It is, 
however, conceivable that general requirements for a versatile and representationally powerful cortical 
circuit architecture are realized by the common design. If this was true, the evolutionary benefit of meeting 
these requirements might have driven the adoption of large-scale self-organization and the emergence of the 
common design over evolutionary times. 

The most prominent candidate for such a general requirement is the hypothesis of a coverage - continuity - 
compromise (e.g. [19,21 55,56]). It states that the columnar organization is shaped to achieve an optimal 
tradeoff between the coverage of the space of visual stimulus features and the continuity of their cortical 
representation. On the one hand, each combination of stimulus features should be well-represented in a 
cortical map to avoid "blindness" to stimuli with particular feature combinations. On the other hand, the 
wiring cost to establish connections within the map of orientation preference should be kept low. This can 
be achieved if neurons that are physically close in the cortex tend to have similar stimulus preferences. 
These two design goals generally compete with each other. The better a cortical representation covers the 
stimulus space, the more discontinuous it has to be. The tradeoff between the two aspects can be modeled 
in what is called a dimension reduction framework in which cortical maps are viewed as two-dimensional 
sheets which fold and twist in a higher-dimensional stimulus space (see Fig. [T]) to cover it as uniformly as 



possible while minimizing some measure of continuity 21 57 58 . 
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From prior work, the coverage continuity compromise appears to be a promising candidate for a principle 
to explain visual cortical functional architecture. Firstly, many studies have reported good qualitative 
agreement between the layout of numerically obtained dimension reducing maps and experimental 



observations 19,21 42,55-57,59 71 . Secondly, geometric relationships between the representations of 
different visual features such as orientation, spatial frequency and ocular dominance have been reproduced 



by dimension reduction models ^25 56,63 65 68,72 73 



Mathematically, the dimension reduction hypothesis implies that the layouts of cortical maps can be 
understood as optima or near optima (global or local minima) of a free energy functional which penalizes 
both "stimulus scotomas" and map discontinuity. Unfortunately, there is currently no dimension reduction 
model for which the precise layouts of optimal or nearly optimal solutions have been analytically 
established. To justify the conclusion that the tradeoff between coverage and continuity favors the common 
rules of 0PM design found in experiment, knowledge of optimal dimension-reducing mappings however 
appears essential. 

Precise knowledge of the spatial organization of optimal and nearly optimal mappings is also critical for 
distinguishing between optimization theories and frozen noise scenarios of visual cortical development. In a 
frozen noise scenario, essentially random factors such as haphazard wiring 74 , the impact of spontaneous 
activity patterns 75 or an idiosyncratic set of visual experiences 76 determine the emerging pattern of 
preferred orientations. This pattern is then assumed to be "frozen" by an unknown mechanism, capable of 
preventing further modification of preferred orientations by ongoing synaptic turnover and 
activity-dependent plasticity. Conceptually, a frozen noise scenario is diametrically opposed to any kind of 
optimization theory. Even if the reorganization of the pattern prior to freezing was to follow a gradient 
descent with respect to some cost function, the early stopping implies that the layout is neither a local nor 
a global minimum of this functional. Importantly, the layout of transient states is known to exhibit 
universal properties that can be completely independent of model details 25','32l. As a consequence, an 
infinite set of distinct optimization principles will generate the same spatial structure of transient states. 
This implies in turn that the frozen transient layout is not specifically shaped by any particular 
optimization principle. Map layout will thus in principle only be informative about design or optimization 
principles of cortical processing architectures if maps are not just frozen transients. 

In practice, however, the predictions of frozen noise and optimization theories might be hard to distinguish. 
Ambiguity between these mutually exclusive theories would result in particular, if the energy landscape of 
the optimization principle is so "rugged" that there is essentially a local energy minimum next to any 
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relevant random arrangement. Dimension reduction models are conceptually related to combinatorial 
optimization problems like the Traveling Salesman Problem and many such problems are believed to 
exhibit rugged energy landscapes [77|j79|. It is therefore essential to clarify, whether paradigmatic 
dimension reduction models are characterized by a rugged or a smooth energy landscape. If their energy 
landscapes were smooth with a small number of well-separated local minima, their predictions would be 
easy to distinguish from those of a frozen noise scenario. 

In this study, we examine the classical example of a dimension reduction model, the elastic network (EN) 
model. Since the seminal work of Durbin & Mitchison [21j, the EN model has been widely used to study 



visual cortical representations 25,42,62-65,69 72 80 . The EN model possesses an explicit energy 
functional which trades off a matching constraint which matches cortical cells to particular stimulus 
features via Hebbian learning, with a continuity constraint that minimizes euclidean differences in feature 
space between neighboring points in the cortex [6^. In two ways, the EN model's explicit variational 
structure is very appealing. Firstly, coverage and continuity appear as separate terms in the free energy 
which facilitates the dissection of their relative influences. Secondly, the free energy allows for the 
formulation of a gradient descent dynamics. The emergence of cortical selectivity patterns and their 
convergence towards a minimal energy state in this dynamics might serve as a model for an optimization 
process taking place in postnatal development. 

Following Durbin and Mitchison, we consider the EN model for the joint mapping of two visual features: 
(i) position in visual space, represented in a retinotopic map (RM) and (ii) line orientation, represented in 
an 0PM. To compute optimal dimensional reducing mappings, we first develop an analytical framework for 
deriving closed-form expressions for fixed points, local minima, and optima of arbitrary optimization 
models for the spatial layout of OPMs and RMs in which predicted maps emerge by a supercritical 
bifurcation as well as for analyzing their stability properties. By applying this framework to different 
instantiations of the EN model, we systematically disentangle the effects of individual model features on 
the repertoire of optimal solutions. We start with the simplest possible model version, a fixed uniform 
retinotopy and an orientation stimulus ensemble with only a single orientation energy and then relax the 
uniform retinotopy assumption incorporating retinotopic distortions. An analysis for a second widely-used 
orientation stimulus ensemble including also unoriented stimuli is given in Appendix I. Surprisingly, in all 
cases, our analysis yields pinwheel-free orientation stripes or stereotypical square arrays of pinwheels as 
local minima or optimal orientation maps of the EN model. Numerical simulations of the EN confirm these 
findings. They indicate, that more complex spatially aperiodic solutions are not dominant and that the 
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energy landscape of the EN model is rather smooth. Our results demonstrate that while aperiodic 
stationary states exist, they are generally unstable in the considered model versions. 
To test whether the EN model is in principle capable of generating complex spatially aperiodic optimal 
orientation maps, we then perform a comprehensive unbiased search of the EN optima for arbitrary 
orientation stimulus distributions. We identify two key parameters determining pattern selection: (i) the 
intracortical interaction range and (ii) the fourth moment of the orientation stimulus distribution function. 
We derive complete phase diagrams summarizing pattern selection in the EN model for fixed as well as 
variable retinotopy. Small interaction ranges together with low to intermediate fourth moment values lead 
to pinwheel-free orientation stripes, rhombic or hexagonal crystalline orientation map layouts as optimal 
states. In the regime of large interaction ranges together with higher fourth moment, we find irregular 
aperiodic 0PM layouts with low pinwheel densities as optima. Only in an extreme and previously 
unconsidered parameter regime of very large interaction ranges and stimulus ensemble distributions with 
an infinite fourth moment, optimal 0PM layouts in the EN model resemble experimentally observed 
aperiodic pinwheel-rich layouts and quantitatively reproduce the recently described species-invariant 
pinwheel statistics. Unexpectedly, we find that the stabilization of such layouts is not achieved by an 
optimal tradeoff between coverage and continuity of a localized population encoding by the maps but 
results from effectively suppressive long-range intracortical interactions in a spatially distributed 
representation of localized stimuli. 

We conclude our re-examination of the EN model with a comparison between different numerical schemes 
for the determination of optimal or nearly optimal mappings. For large numbers of stimuli, numerically 
determined solutions match our analytical predictions, irrespectively of the computational method used. 

Results and Discussion 

Model definition and model symmetries 

We analyze the elastic net (EN) model for the joint optimization of position and orientation selectivity as 
originally introduced by Durbin & Mitchison 21 . In this model, the retinotopic map (RM) is represented 
by a mapping R(x) = (i?i(x), i?2(x)) which describes the receptive field center position of a neuron at 
cortical position x. Any retinotopic map can be decomposed into an affine transformation x X from 
cortical to visual field coordinates, on which a vector-field of retinotopic distortions r(x) is superimposed, 
i.e.: 

R(x) = X + r(x) 
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Figure 2: The elastic net (EN) model, (a) Example orientation preference map (0PM) (color code) 
together with a uniform map of visual space (RM) (grid lines) (b) Position = (sx^Sy) and orientation 
6> of a 'pointlike' stimulus, (c) Cortical activity, evoked by the stimulus in b for the model maps in a. 
Dark regions are activated. Note, that in contrast to SOFM models, the activity pattern does not exhibit 
a stereotypical Gaussian shape, (d) Modification of orientation preference and retinotopy, caused by the 
stimulus in b. Orientation preferences prior to stimulus presentation are indicated with grey bars, after 
stimulus presentation with black bars. Most strongly modified preferences correspond to thick black bars. 
Modifications of orientation preferences and retinotopy are displayed on an exaggerated scale for illustration 
purposes. 
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with appropriately chosen units for x and R. 

The orientation preference map (0PM) is represented by a second complex- valued scalar field ^(x). The 
pattern of orientation preferences 'i^(x) is encoded by the phase of z(x) via 

^(x) = i arg(2:(x)) . 

The absolute value |^(x)| is a measure of the average cortical selectivity at position x. Solving the EN 
model requires to find pairs of maps {r(x),z(x)} that represent an optimal compromise between stimulus 
coverage and map continuity. This is achieved by minimizing a free energy functional 

T^G^C^n (1) 

in which the functional C measures the coverage of a stimulus space and the functional IZ the continuity of 
its cortical representation. The stimulus space is defined by an ensemble {S} of idealized pointlike stimuli, 
each described by two features: = \sz\e^'^^ and = {sx^^Sy) which specify the orientation 6 of the 
stimulus and its position in visual space (Fig. [2b). C and IZ are given by 



= ^/^'^^l|V^(y)||Vr?.^||Vr,(y)||\ 



with V = (dx^dy)^ ^ and rj G [0, 1]. The ratios and cf'^ /rjr control the relative strength of the coverage 
term versus the continuity term for 0PM and RM, respectively. (•••)§ denotes the average over the 
ensemble of stimuli. 

Minima of the energy functional T are stable fixed points of the gradient descent dynamics 

d,z{^) = -2^gljl^F-[z,r](x) (2) 
a.r(x) = -^^.F'-[z,r](x) 

called the EN dynamics in the following. These dynamics read 

dtz{x.) = ([5^ -z(x)]e(x,S,z,r))g +7^ A2;(x) (3) 
S,r(x) = ([s,-X-r(x)]e(x,S,z,r))s+7^, Ar(x), (4) 

where e(x, S, r) is the activity-pattern, evoked by a stimulus S = (s^, Sz) in a model cortex with 
retinotopic distortions r(x) and 0PM z{-k). It is given by 

e(x,...) 



J ^2y g-(|s.-X-r(y)|2)/2a2g-(|s,-^(y)|2)/2a2 ' 



Figure [2] illustrates the general features of the EN dynamics using the example of a single stimulus. Fig. [2^ 
shows a model orientation map with a superimposed uniform representation of visual space. A single 
pointlike, oriented stimulus S = (s^,^^) with position = (sx^Sy) and orientation 6 = l/2arg(s;^) (Fig. 
^) evokes a cortical activity pattern e(x, S, r) (Fig. The activity-pattern in this example is of 
roughly Gaussian shape and is centered at the point, where the location of the stimulus is represented in 
cortical space. However, depending on the model parameters and the stimulus, the cortical activity pattern 
may assume as well a more complex form (see also Discussion). According to Eqs. (|3j|4|, each stimulus 
and the evoked activity pattern induce a modification of the orientation map and retinotopic map, shown 
in Fig. [2]i. Orientation preference in the activated regions is shifted towards the orientation of the 
stimulus. The representation of visual space in the activated regions is locally contracted towards the 
position of the stimulus. Modifications of cortical selectivities occur due to randomly chosen stimuli and 
are set proportional to a very small learning rate. Substantial changes of cortical representations occur 
slowly through the cumulative effect of a large number of activity patterns and stimuli. These effective 
changes are described by the two deterministic equations for the rearrangement of cortical selectivities Eqs. 
(|3j [4| which are obtained by stimulus- averaging the modifications due to single activity patterns in the 
discrete stimulus model 25 . One thus expects that the optimal selectivity patterns and also the way in 
which cortical selectivities change over time are determined by the statistical properties of the stimulus 
ensemble. In the following, we assume that the stimulus ensemble satisfies three properties: (i) The 
stimulus locations are uniformly distributed across visual space, (ii) For the distribution of stimulus 
orientations, l^^l and are independent, (iii) Orientations are distributed uniformly in [0,7r]. 
These conditions are fulfilled by stimulus ensembles used in virtually all prior studies of dimension 



reduction models for visual cortical architecture (e.g. 19}[2H[25 64,65 71,72 80-82]). They imply several 



symmetries of the model dynamics Eqs. ([3j|4]). Due to the first property, the EN dynamics are equivariant 
under translations 

Tyz(x) = z(x + y) 
Tyr(x) = r(x + y), 

rotations 

^/3r(x) = r(l]_/3x) 
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with 2x2 rotation matrix 



and reflections 



cos j3 — sin P 
sin j3 cos j3 



Pz{x.) = z(^x) 
Pr(x) = ^r(^x) , 

where ^ = diag(— 1, 1) is the 2x2 reflection matrix. Equivariance means that 

fyF^[z,r] = F^[fyz,fyr] 
RpF'[z,r] = F^[Rpz,Rpv] 
PF^[2,r] = F'[Pz,Py], 



(5) 
(6) 

(7) 



with mutatis mutandis the same equations fulfilled by the vector-field F^[2:,r]. 

As a consequence, patterns that can be converted into one another by translation, rotation or reflection of 
the cortical layers represent equivalent solutions of the model Eqs. ([3j[4|, by construction. Due to the 
second assumption, the dynamics is also equivariant with respect to shifts in orientation S'0z(x) = e*'^z(x), 
i.e. 



(8) 
(9) 



Thus, two patterns are also equivalent solutions of the model, if their layout of orientation domains and 
retinotopic distortions is identical, but the preferred orientations differ everywhere by the same constant 
angle. 

Without loss of generality, we normalize the ensembles of orientation stimuli such that 

(l"^^l^)s ~ (I'^^P) ~ ^ throughout this paper. This normalization can always be restored by a rescaling of 



z(x) (see 25,69^). 



Our formulation of the dimension-reduction problem in the EN model utilizes a continuum description, 
both for cortical space and the set of visual stimuli. This facilitates mathematical treatment and appears 
appropriate, given the high number of cortical neurons under one square millimeter of cortical surface (e.g. 
roughly 70000 in cat VI 83 ). Even an hypothesized neuronal mono-layer would consist of more than 
20x20 neurons per hypercolumn area A^, constituting a quite dense sampling of the spatial periodicity. 
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Figure 3: The linearization of the EN model dynamics around the unselective fixed point, (a) 

Eigenvalue spectra of the linearized retinotopy dynamics for longitudinal mode (A£(/c), blue trace) and 
transversal mode (A^(/c), red trace), (b) Longitudinal mode ~ k^e*^'^^ + c.c. (left) and transversal mode 
~ k0+7^/2e*^^^ + c.c. (right), (c) Spectrum of eigenvalues of the linearized 0PM dynamics (red trace) for 
a < (J*{r]). Orange region marks the unstable annulus of Fourier modes (critical circle), (d) Stability regions 
of the nonselective state in the EN model. The stability border is given by cr^{r]) (Eq. (Il2|). 



Treating the feature space as a continuum implements the concept that the cortical representation has to 
cover as good as possible the infinite multiplicity of conceivable stimulus feature combinations. 



The orientation unselective fixed point 

Two stationary solutions of the model can be established from symmetry. The simplest of these is the 
orientation unselective state with z(x) =0 and uniform mapping of visual space r(x) = 0. Firstly by the 
shift symmetry Eq. ([8|, we find that z(x) = is a fixed point of Eq. ([3|. Secondly by reflect ional and 
rotational symmetry Eqs. ([5j[7|), we see that the right-hand side of Eq. Q has to vanish and hence the 
orientation unselective state with uniform mapping of visual space is a fixed point of Eqs. (|3j[4|. 
This homogeneous unselective state thus minimizes the EN energy functional if it is a stable solution of 
Eqs. (|3|[4|. The stability can be determined by considering the linearized dynamics of small deviations 
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{r(x),z(x)} around this state. These hnearized dynamics read 

atr(x) ~ Lr[r] = J d^ye-^^Ar{y)+7]rAr{x) (10) 

dtz{x) ^ L,[z]=(^^-l^z{^)+vAz{x)-^ Jd^ye-^ziy). (11) 

where (A)^j = {xi — yi){xj — yj) — 2(j'^Sij with Sij being Kronecker's delta. We first note that the hnearized 
dynamics of retinotopic distortions and orientation preference decouple. Thus, up to linear order and near 
the homogeneous fixed point, both cortical representation evolve independently and the stability properties 
of the unselective state can be obtained by a separate examination of the stability properties of both 
cortical representations. 

The eigenfunctions of the linearized retinotopy dynamics L^[r] can be calculated by Fourier-transforming 



Eq. ([10|): 

2 

where /c = |k| and z = 1,2. A diagonalization of this matrix equation yields the eigenvalues 

= -k\r]r + e-^'^V^), A^. = -rjrk^ 

with corresponding eigenfunctions (in real space) 

rL(x) = k^e*^^"" + c.c. 
rT(x) = k^+^/se^^^^ + c.c. , 

where k^ = |k|(cos (/), sin (/))-^. These eigenfunctions are longitudinal (L) or transversal (T) wave patterns. 
In the longitudinal wave, the retinotopic distortion vector r(x) lies parallel to k which leads to a 
"compression wave" (Fig. [sJd, left). In the transversal wave pattern (Fig. [sJd, right), the retinotopic 
distortion vector is orthogonal to k. We note that the linearized Kohonen model [6l] was previously found 
to exhibit the same set of eigenfunctions 81 . Because both spectra of eigenvalues A^, A£ are smaller than 
zero for every a > 0, 77^ > 0, and k > (Fig. [S^i), the uniform retinotopy r(x) = is a stable fixed point of 
Eq. Q irrespective of parameter choice. 

The eigenfunctions of the linearized 0PM dynamics [z] are Fourier modes ~ e*^^ by translational 
symmetry. By rotational symmetry, their eigenvalues only depend on the wave number k and are given by 



X'{k) = -1 + (1 - e-^'-') - r]k^ . 
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(see 25 ). This spectrum of eigenvalues is depicted in Fig. For 7? > 0, X^{k) has a single maximum at 
fee = ^^/WW- For 



a > a*{r]) ^/l^r]\nr]- rj (12) 

this maximal eigenvalue r = X^ikc) is negative. Hence, the unselective state with uniform retinotopy is a 
stable fixed point of Eqs. ([sj |4| and the only known solution of the EN model in this parameter range. 
For cr < cr* (77), the maximal eigenvalue r is positive, and the nonselective state is unstable with respect to a 
band of Fourier modes ~ e*^^ with wave numbers around |k| ^ kc (see Fig. [s]^). This annulus of unstable 
Fourier modes is called the critical circle. The finite wavelength instability |84||86] (or Turing 
instability [87]) leads to the emergence of a pattern of orientation preference with characteristic spacing 
A = 2i\ jkc from the nonselective state on a characteristic timescale r = 1/r. 



One should note that as in other models for the self-organization of orientation columns, e.g. 15 57 , the 
characteristic spatial scale A arises from effective intracortical interactions of 'Mexican-hat' structure 
(short-range facilitation, longer-ranged suppression). The short-range facilitation in the linearized EN 



dynamics is represented by the first two terms on the right hand side of Eq. (11). Since a < 1 in the 
pattern forming regime, the prefactor in front of the first term is positive. Due to the second, Laplacian 
term, it is favored that neighboring units share selectivity properties, a process mediated by short-range 



facilitation. Longer-ranged suppression is represented by the convolution term in Eq. (11). 
Mathematically, this term directly results from the soft-competition in the "activity-dependent" coverage 
term of Eq. ([T]). The local facilitation is jointly mediated by coverage (first term) and continuity (second 
term) contributions. 

Fig. [sji summarizes the result of the linear stability analysis of the nonselective state. For a > cF*(rf)^ the 
orientation unselective state with uniform retinotopy is a minimum of the EN free energy and also the 
global minimum. For cr < cr*{r]) this state represents a maximum of the energy functional and the minima 
must thus exhibit a space-dependent pattern of orientation selectivities. 



Orientation stripes 

Within the potentially infinite set of orientation selective fixed points of the model, one class of solutions 
can be established from symmetry: {r(x) = 0, z(x) = A^e^^^}. In these pinwheel-free states, orientation 
preference is constant along one axis in cortex (perpendicular to the vector k), and each orientation is 
represented in equal proportion (see Fig. |4^). Retinotopy is perfectly uniform. Although this state may 
appear too simple to be biologically relevant, we will see that it plays a fundamental role in the state space 
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Figure 4: Exact and approximate orientation selective fixed points of OPM optimization models 

(a) Pinwheel free orientation stripe (OS) pattern. Diagram shows the position of the wave vector in Fourier 
space, (b) Rhombic Pinwheel crystal (rPWC) with four nonzero wave vectors, (c) Essentially complex 
planforms (ECPs). The index n indicates the number of nonzero wave vectors. The index i enumerates 
nonequivalent configurations of wave vectors with the same n, starting with i = for the most anisotropic 
planform. For n = 3, 5, and 15, there are 2, 4, and 612 different ECPs, respectively. OPM layouts become 
more irregular with increasing n. 
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of the EN model. It is therefore useful to establish its existence and basic characteristics. The existence of 
orientation stripe (OS) solutions follows directly from the model's symmetries Eqs. (5][9). Computing 



demonstrates that F^[e*^^,0] is proportional to e*^^. This establishes that the subspace of functions 
~ g*kx invariant under the dynamics Eq. (|3|. For Aq = 0, we recover the trivial fixed point of the EN 
dynamics by construction, as shown above. This means that within this subspace = is either a 
minimum or a maximum of the EN energy functional Eq. ([T]). Furthermore, for Aq ^ oo the EN energy 
tends to infinity. If the trivial fixed point is unstable, it corresponds to a maximum of the EN energy 
functional. Therefore, there must exist at least one minimum with Aq ^ in the subspace of functions 
~ e*^^ which then corresponds to a stationary state of the EN dynamics. 

Regarding the dynamics of retinotopic deviations, the model's symmetries Eqs. ([5][9| can be invoked to 
show that for the state {O^Aqc^^^}^ the right-hand side of Eq. Q has to be constant in space: 

Ty[F^[e^^^,0]] = F^[Ty[e^^^],Ty[0]] = F"^ [e'^^ e'^"" , 0] = F^[e^^^,0] 

If this constant was non-zero the retinotopic map would drift with constant velocity. This, however, is 
impossible in a variational dynamics such that this constant must vanish. The orientation stripe solution 
(Fig. [4|l) is to the best of our knowledge the only exact nontrivial stationary solution of Eqs. (|3j [4| that 
can be established without any approximations. 



Doubly periodic and quasiperiodic solutions 

In the EN model as considered in this study, the maps of visual space and orientation preference are jointly 
optimized to trade off coverage and continuity leading to mutual interactions between the two cortical 
representations. These mutual interactions vanish in the rigid retinotopy limit 77^^00 and the perfectly 
uniform retinotopy becomes an optimal solution for arbitrary orientation column layout z{'x.). As it is not 
clear how essential the mutual interactions with position specificity are in shaping the optimal orientation 
column layout, we continue our investigation of solution classes by considering global minima of 
optimization models with fixed uniform retinotopy. The mutual interactions will be taken into account in a 
subsequent step. 

In the rigid retinotopy limit, minima of the energy functional are stable stationary states of the dynamics 
of the orientation preference map Eq. ([3| with r(x) = 0. To compute orientation selective stationary 



16 



solutions of this 0PM dynamics, we employ that in the vicinity of a supercritical bifurcation where the 
nonselective fixed point z(x) = becomes unstable, the entire set of nontrivial fixed points is determined 
by the third-order terms of the Volterra series representation of the operator 0] [35l|85[|86[ |88|. The 
symmetries Eqs. ([5][9| restrict the general form of such a third-order approximation for any model of 0PM 
optimization to 

dtz{^)^L,[z]+Ni[z,z,z], (13) 



where the cubic operator is written in trilinear form, i.e. 



7V| 



j k I 



i,Zk, zi\ . 



In particular, all even terms in the Volterra Series representation of F^[z^ 0] vanish due to the 

Shift- Symmetry Eqs. ([8j|9|. Explicit analytic computation of the cubic nonlinearities for the EN model is 

cumbersome but not difficult (see Methods) and yields a sum 

11 

7V|[z,z,z] = ^a,7V^'[z,z,z]. (14) 

i=i 

The individual nonlinear operators are with one exception nonlocal convolution-type operators and are 



given in the Methods section (Eq. (38)), together with a detailed description of their derivation. Only the 



coefficients aj depend on the properties of the ensemble of oriented stimuli. 



To calculate the ffxed points of Eq. (13), we use a perturbative method called weakly nonlinear analysis 



that enables us to analytically examine the structure and stability of inhomogeneous stationary solutions in 



the vicinity of a ffnite- wavelength instability. Here, we examine the stability of so-called planforms 84 - 86 
Planforms are patterns that are composed of a finite number of Fourier components, such as 



^(x) = ^A,(t)e''^'- 

j 



for a pattern of orientation columns. With the above planform ansatz, we neglect any spatial dependency 
of the amplitudes Aj (t) for example due to long- wave deformations for the sake of simplicity and analytical 
tractability. When the dynamics is close to a finite wavelength instability, the essential Fourier components 
of the emerging pattern are located on the critical circle |kj| = kc. The dynamic equations for the 
amplitudes of these Fourier components are called amplitude equations. For a discrete number of N 
Fourier components of z{x.) whose wave vectors lie equally spaced on the critical circle, the most general 
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system of amplitude equations compatible with the symmetries Eqs. ([5j|9|) has the form [35 

N N 

Ai = rAi — Ai ^ ^ Qij I Aj I ^ — A^- ^ ^ fijAjAj- , (15) 

with r > 0. Here, gij and /^^ are the real- valued coupling coefficients between the amplitudes Ai and Aj. 
They depend on the differences between indices \i — j\ and are entirely determined by the nonlinearity 



N^[z^z^z] in Eq. (13). If the wave vectors = (cosa^, sina^j/cc are parametrized by the angles a^, then 
the coefficients gij and fij are functions only of the angle a = \ai — aj \ between the wave vectors and 
kj. One can thus obtain the coupling coefficients from two continuous functions g{a) and f{a) that can be 
obtained from the nonlinearity N^[z,z^z] (see Methods for details). In the following, these functions are 



called angle-dependent interaction functions. The amplitude equations Eq. ( 15 ) are variational if and only 
if gij and fij are real-valued. In this case they can be derived via 

from an energy 

N 1 ^ 1 ^ 

Ua = -rJ2\A,f + 9^J\A.\'\AJf + 2 E h^A-AjA^- . (16) 

i=l i,j=l i,j=l 

If the coefficients gij and fij are derived from Eq. ([T]), the energy Ua for a given planform solution 
corresponds to the energy density of the EN energy functional considering only terms up to fourth order in 
z(x). 



The amplitude equations Eq. ( 15 ) enable to calculate an infinite set of orientation selective fixed points. 
For the above orientation stripe solution with one nonzero wave vector z{'k) = Aqg*^^, the amplitude 
equations predict the so far undetermined amplitude 

l^ol' = f (17) 

and its energy 

Since ga > 0, this shows that OS stationary solutions only exist for r > 0, i.e. in the symmetry breaking 
regime. As for all following fixed-points, Uqs specifies the energy difference to the homogeneous unselective 
state z(x) = 0. 

A second class of stationary solutions can be found with the ansatz 

z(x) = Aie'^'"" + ^26^^^^ + Ase-'^'"" + A^e-'^'"" , 
18 



with amplitudes Aj = \ Aj\e^^^ and Z(ki,k2) = a > 0. By inserting this ansatz into Eq. (15) and assuming 
uniform amphtude \Ai \ = \A2\ = \A2\ = {A^l = we obtain 

= ^ • (19) 

The phase relations of the four amplitudes are given by 

01+03 00 
02+04 = 00 + TT . 

These solutions describe a regular rhombic lattice of pinwheels and are therefore called rhombic pinwheel 
crystals (rPWC) in the following. Three phases can be chosen arbitrarily according to the two above 
conditions, e.g. 0o, Aq = 0i — 03 and Ai = ^2 — 04- For an rPWC parametrized by these phases, Aq shifts 
the absolute positions of the pinwheels in x-direction, Ai shifts the absolute positions of the pinwheels in 
y-direction, and 0o shifts all the preferred orientations by a constant angle. The energy of an rPWC 
solution is 

2r 

^rPWC = ^ . (20) 

^00 + gOir + gOa + gOn-a " ^JOa 

An example of such a solution is depicted in Fig. [4b. We note that rPWCs have been previously found in 



several other models for 0PM development 27,31 37,39,89 . The pinwheel density p of an rPWC, i.e. the 
number of pinwheels in an area of size A^, is equal to p = 4 sin a 54 . The angle a which minimizes the 



energy Urpwc can be computed by maximizing the function 

s{a) = goa + goTT-a - 2/oa (21) 



in the denominator of Eq. ( 20 ) 



The two solution classes discussed so far, namely OS and rPWCs, exhibit one prominent feature, absent in 
experimentally observed cortical OPMs, namely perfect spatial periodicity. Many cortical maps including 
OPMs do not resemble a crystal-like grid of repeating units. Rather the maps are characterized by roughly 
repetitive but aperiodic spatial arrangement of feature preferences (e.g. [sjjlO]). This does not imply that 
the precise layout of columns is arbitrary. It rather means that the rules of column design cannot be 
exhaustively characterized by mapping a "representative" hypercolumn. 

Previous studies of abstract models of 0PM development introduced the family of so-called essentially 



complex planforms (ECPs) as stationary solutions of Eq. (15). This solution class encompasses a large 



variety of realistic quasiperiodic 0PM layouts and is therefore a good candidate solution class for models of 
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0PM layout. In addition, Kaschube et al. demonstrated that models in which these are optimal solutions 
can reproduce all essential features of the common 0PM design in ferret, tree-shrew, and galago |38j . 
An n-ECP solution can be written as 

n 

with n = N/2 wave vectors kj = /Cc(cos(7rj/n), sin(7rj/n)) distributed equidistantly on the upper half of the 
critical circle, complex amplitudes Aj and binary variables Ij = ±1 determining whether the mode with 
wave vector kj or — kj is active (nonzero). Because these planforms cannot realize a real- valued function 



they are called essentially complex [35]. For an n-ECP, the third term on the right hand side of Eq. (15) 
vanishes and the amplitude equations for the active modes Ai reduce to a system of Landau equations 

n 

A-i = vAi — Ai ^ ^ 9ij I Aj I . 
where gij is the n x n-coupling matrix for the active modes. Consequently, the stationary amplitudes obey 

n 

|A,p = r^(g-^),,. (22) 

The energy of an n-ECP is given by 

t^ECP = -^E(s"%- (23) 
We note that this energy in general depends on the configuration of active modes, given by the /j's, and 
therefore planforms with the same number of active modes may not be energetically degenerate. 
Families of n-ECP solutions are depicted in Fig. [4J3. The 1-ECP corresponds to the pinwheel-free OS 
pattern discussed above. For fixed n > 3, there are multiple planforms not related by symmetry operations 
which considerably differ in their spatial layouts. For n > 4, the patterns are spatially quasiperiodic, and 
are a generalization of the so-called Newell- Pomeau turbulent crystal 90,91]. For n > 10, their layouts 
resemble experimentally observed OPMs. Different n-ECPs however differ considerably in their pinwheel 
density. Planforms whose nonzero wave vectors are distributed isotropically on the critical circle typically 
have a high pinwheel density (see Fig. [4]^, n=15 lower right). Anisotropic planforms generally contain 
considerably fewer pinwheels (see Fig. |4]3, n=15 lower left). All large n-ECPs, however, exhibit a complex 
quasiperiodic spatial layout and a nonzero density of pinwheels. 

In order to demonstrate that a certain planform is an optimal solution of an optimization model for 0PM 
layouts in which patterns emerge via a supercritical bifurcation, we not only have to show that it is a 
stationary solution of the amplitude equations but have to analyze its stability properties with respect to 
the gradient descent dynamics as well as its energy compared to all other candidate solutions. 
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Many stability properties can be characterized by examining the amphtude equations (15). In principle, 
the stability range of an n-ECPs may be bounded by two different instability mechanisms: (i) an intrinsic 
instability by which stationary solutions with n active modes decay into ones with lower n. (ii) an extrinsic 
instability by which stationary solutions with a "too low" number of modes are unstable to the growth of 
additional active modes. These instabilities can constrain the range of stable n to a small finite set around 



a typical n 35 88 . A mathematical evaluation of both criteria leads to precise conditions for extrinsic and 
intrinsic stability of a planform (Methods). In the following, a planform is said to be stable, if it is both 
extrinsically and intrinsically stable. A planform is said to be an optimum (or optimal solution) if it is 
stable and possesses the minimal energy among all other stationary planform solutions. 
Taken together, this amplitude equation approach enables to analytically compute the fixed points and 
optima of arbitrary optimization models for visual cortical map layout in which the functional architecture 
is completely specified by the pattern of orientation columns z(x) and emerges via a supercritical 
bifurcation. Via a third-order expansion of the energy functional together with weakly nonlinear analysis, 
the otherwise analytically intractable partial integro-differential equation for 0PM layouts reduces to a 
much simpler system of ordinary differential equations, the amplitude equations. Using these, several 
families of solutions, orientation stripes, rhombic pinwheel crystals and essentially complex planforms, can 
be systematically evaluated and comprehensively compared to identify sets of unstable, stable and optimal, 
i.e. lowest energy fixed points. 

As already mentioned, the above approach is suitable for arbitrary optimization models for visual cortical 
map layout in which the functional architecture is completely specified by the pattern of orientation 
columns z(x) which in the EN model is fulfilled in the rigid retinotopy limit. We now start by considering 
the EN optimal solutions in this limit and subsequently generalize this approach to models in which the 
visual cortical architecture is jointly specified by maps of orientation and position preference that are 
matched to one another. 

Representing an ensemble of "bar" -stimuli 

We start our investigation of optimal dimension-reducing mappings in the EN model using the simplest 
and most frequently used orientation stimulus ensemble, the distribution with 5;^-values uniformly arranged 
on a rme; with radius r,^ = ^2 [57l[64}|66)|92] . We call this stimulus ensemble the circular stimulus 
ensemble in the following. According to the linear stability analysis of the nonselective fixed point, at the 



point of instability, we choose a = cr*{r]) such that the linearization Eq. (11) is completely characterized by 
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Figure 5: Angle-dependent interaction functions for the EN model with fixed retinotopy and 
circular orientation stimulus ensemble. (a,b) g{a) and f{a) for a/A = 0.1 (a) and a/A = 0.2 (b). 



the continuity parameter rj. Equivalent to specifying 77 is to fix the ratio of activation range a and column 
spacing A 

a/A=^y/]^i{lM (24) 

as a more intuitive parameter. This ratio measures the effective interaction-range relative to the expected 
spacing of the orientation preference pattern. In abstract optimization models for 0PM development a 



similar quantity has been demonstrated to be a crucial determinant of pattern selection 35 , 88 . We note 



however, that due to the logarithmic dependence of a/A on 77, a slight variation of the effective interaction 
range may correspond to a variation of the continuity parameter rj over several orders of magnitude. 
In order to investigate the stability of stationary planform solutions in the EN model with a circular 
orientation stimulus ensemble, we have to determine the angle-dependent interaction functions g{a) and 



f{a). For the coefficients aj in Eq. (14) we obtain 



ai = 
a4 = 
ay = 
aio 
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The angle-dependent interaction functions of the EN model with a circular orientation stimulus ensemble 
are then given by 



9{(^) = 



2e '"^^ — e 



2/c^cr^(cos 



a-1) ^ 



sinh'^ {l/2kl(j^ cos a) 



-kta' 



cr'^ cos a 



(25) 



(cosh(2/c^a^ cos a) + 2 cosh(/c^a^ cos a)) + 2e 
+ ^ (e-^^'^'cosh(2A:^a^cosa) - l) + ^e'^ 

These functions are depicted in Fig. [sjfor two different values of the interaction range a/A. We note that 
both functions are positive for all a/ A which is a sufficient condition for a supercritical bifurcation from 
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the homogeneous nonselective state in the EN model. 



Finally, by minimizing the function s{a) in Eq. (21), we find that the angle a which minimizes the energy 



of rPWC fixed-point is a = 7r/2. This corresponds to a square array of pinwheels (sPWC). Due to the 
orthogonal arrangement oblique and cardinal orientation columns and the maximized pinwheel density of 
p = 4, the square array of pinwheels has the maximal coverage among all rPWC solutions. 

Optimal solutions close to the pattern formation threshold 

We first tested for the stability of pinwheel-free OS solutions and the sPWCs, by analytical evaluation of 
the criteria for intrinsic and extrinsic stability (Methods). We found both^ OS and sPWCs, to be 
intrinsically and extrinsically stable for all a/ A. Next, we tested for the stability of n-ECP solutions with 
2 < n < 20. We found all n-ECP configurations with 2 < n < 20 to be intrinsically unstable for all <j/A. 
Hence, none of these planforms represent optimal solutions of the EN model with a circular stimulus 
ensemble, while both OS and sPWC are always local minima of the energy functional. 



By evaluating the energy assigned to the sPWC (Eq. ( |20| )) and OS (Eq. (18)), we next identified two 
different regimes: (i) For short interaction range a/A < 0.122 the sPWC possesses minimal energy and is 
therefore the predicted global minimum, (ii) For a/A > 0.122 the OS is optimal. Figure [6|i shows the 
resulting simple phase diagram. The sPWC and OS are separated by a phase border at a/A ?^ 0.122. We 
numerically confirmed these analytical predictions by extensive simulations of Eq. (|3| with r(x) = and 
the circular stimulus ensemble (see Methods for details). Fig. [gJd-c show snapshots of a representative 
simulation with short interaction range (r = 0.1, a/A = 0.1 (?7 = 0.67)). After the phase of initial pattern 
emergence (symmetry breaking), the 0PM layout rapidly approaches a square array of pinwheels, the 
analytically predicted optimum (Fig. |6]i). Pinwheel density time courses (see Methods) display a rapid 
convergence to a value close to the predicted density of 4 (Fig. [6^). Fig. |6]F shows the stationary mean 
squared amplitudes of the pattern obtained for different values of the control parameter r (black circles). 



For small control parameters, the pattern amplitude is perfectly predicted by Eq. (19) (solid green line). 
Fig. |6^-h show snapshots of a typical simulation with longer interaction range (r = 0.1, 
a/A = 0.15 (?7 = 0.41)). After the emergence of an 0PM with numerous pinwheels, pinwheels undergo 
pairwise annihilation as previously described for various models of 0PM development and 



optimization 25 27,35 . The OP pattern converges to a pinwheel-free stripe pattern, which is the 
analytically computed optimal solution in this parameter regime (Fig. [gJ). Pinwheel densities decay 
towards zero over the time course of the simulations (Fig. |6]). Also in this parameter regime, the mean 
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Figure 6: Optimal solutions of the EN model with a circular orientation stimulus ensemble 
|57y64| - [66H 92 1 and fixed representation of visual space, (a) At criticality, the phase space of this model 
is parameterized by either the continuity parameter r] (blue labels) or the interaction range a/A (red labels, 
see text), (b-c) OPMs (b) and their power spectra (c) in a simulation of Eq. (|3| with r(x) = and r = 0.1, 
cr/A = 0.1 (t^ = 0.67) and circular stimulus ensemble, (d) Analytically predicted optimum for a/A < 0.122 
(quadratic pinwheel crystal), (e) Pinwheel density time courses for four different simulations (parameters 
as in b; gray traces, individual realizations; black trace, simulation in b; red trace, mean value), (f) Mean 
squared amplitude of the stationary pattern, obtained in simulations (parameters as in b) for different values 
of the control parameter r (black circles) and analytically predicted value (solid green line), (g-h) OPMs 
(g) and their power spectra (h) in a simulation of Eq. ([3| with r(x) = and a/A = 0.15 (t^ = 0.41) (other 
parameters as in b). (i) Analytically predicted optimum for a/A > 0.122 (orientation stripes), (j) Pinwheel 
density time courses for four different simulations (parameters as in g; gray traces, individual realizations; 
black trace, simulation in g; red trace, mean value), (k) Mean squared amplitude of the stationary pattern, 
obtained in simulations (parameters as in g) for different values of the control parameter r (black circles) 
and analytically predicted value (solid green line). 
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squared amplitude of the pattern is well-predicted Eq. ( |17[ ) for small r. (Fig. [6}^). 

In summary, the phase diagram of the EN model with a circular stimulus ensemble close to threshold is 
divided into two regions: (i) for a small interaction range (large continuity parameter) a square array of 
pinwheels is the optimal dimension-reducing mapping and (ii) for a larger interaction range (small 
continuity parameter) orientation stripes are the optimal dimension-reducing mapping. Both states are 
stable throughout the entire parameter range. All other planforms, in particular quasiperiodic n-ECPs are 
unstable. 

At first sight, this structure of the EN phase diagram may appear rather counterintuitive. A solution with 
many pinwheel-defects is energetically favored over a solution with no defects in a regime with large 
continuity parameter where discontinuity should be strongly penalized in the EN energy term. However, a 
large continuity parameter at pattern formation threshold inevitably leads to a short interaction range a 
compared to the characteristic spacing A (see Eq. ([24|). In such a regime, the gain in coverage by 
representing many orientation stimuli in a small area spanning the typical interaction range, e.g. with a 
pinwheel, is very high. Our results show that the gain in coverage by a spatially regular positioning of 
pinwheels outweighs the accompanied loss in continuity above a certain value of the continuity parameter. 

EN dynamics far from pattern formation threshold 

Close to pattern formation threshold, we found only two stable solutions, namely OS and sPWCs. Neither 
of the two exhibits the characteristic aperiodic and pinwheel-rich organization of experimentally observed 
orientation preference maps. Furthermore, the pinwheel densities of both solutions (p = for OS and p = 4 



for sPWCs) differ considerably from experimentally observed values 38 around 3.14. One way towards 
more realistic stable stationary states might be to increase the distance from pattern formation threshold. 
In fact, further away from threshold, our perturbative calculations may fail to correctly predict optimal 
solutions of the model due to the increasing influence of higher order terms in the Volterra series expansion 
of the right hand side in Eq. 

To asses this possibility, we simulated Eq. Q with r(x) = and a circular stimulus ensemble for very large 
values of the control parameter r. Fig. [7] displays snapshots of such a simulation for r = 0.8 as well as their 
pinwheel density time courses for two different values of <j/A. Pinwheel annihilation in the case of large 
cr/ A is less rapid than close to threshold (Fig. [7^, b). The 0PM nevertheless converges towards a layout 
with rather low pinwheel density with pinwheel- free stripe-like domains of different directions joined by 
domains with essentially rhombic crystalline pinwheel arrangement. The linear zones increase their size 
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Figure 7: Numerical analysis of the EN dynamics with circular orientation stimulus ensemble 
and fixed representation of visual space far from pattern formation threshold, (a) OPMs and 
their power spectra in a simulation of Eq. (|3| with r(x) = 0, r = 0.8, a/ A = 0.3 (r^ = 0.028) and circular 
orientation stimulus ensemble. Pinwheel density time courses for four different simulations (parameters as 
in a; gray traces, individual realizations; black trace, simulation in a; red trace, mean value) (c, d) OPMs 
and their power spectra in a simulation of Eq. ([3| with r(x) = 0, r = 0.8, a/ A = 0.12 (t^ = 0.57) and circular 
orientation stimulus ensemble, (d) Pinwheel density time courses for four different simulations (parameters 
as in c; gray traces, individual realizations; black trace, simulation in c; red trace, mean value) 
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over the time course of the simulations, eventuahy leading to stripe-patterns for large simulation times. For 
smaller interaction ranges a/ A, the 0PM layout rapidly converges towards a crystal-like rhombic 
arrangement of pinwheels, however containing several dislocations (Fig. [24^ ) [85]. Dislocations are defects 
of roll or square patterns, where two rolls or squares merge into one, thus increasing the local wavelength of 



the pattern 84,86 . Nevertheless, for all simulations, the pinwheel density rapidly reaches a value close to 4 
(Fig. ^) and the square arrangement of pinwheels is readily recognizable. Both features, the dislocations 
in the rhombic patterns and domain walls in the stripe patterns, have been frequently observed in 



pattern- forming systems far from threshold 85 , 86 



In summary, the behavior of the EN dynamics with circular stimulus ensemble far from pattern formation 
threshold agrees very well with our analytical predictions close to threshold. Again, orientation stripes and 
quadratic pinwheel crystals are identified as the only stationary solutions. Aperiodic and pinwheel-rich 
patterns which resemble experimentally 0PM layouts were not observed. 

Taking retinotopic distortions into account 

So far, we have examined the optimal solutions of the EN model for the simplest and most widely used 
orientation stimulus ensemble. Somewhat unexpected from previous reports, the optimal states in this case 
do not exhibit the irregular structure of experimentally observed orientation maps. Our treatment however 
differs from previous approaches in that the mapping of visual space so far was assumed to be undistorted 
and fixed, i.e. r(x) = 0. We recall that in their seminal publication, Durbin and Mitchison in particular 
demonstrated interesting correlations between the map of orientation preference and the map of visual 
space 



21 . These correlations suggest a strong coupling between the two that may completely alter the 
model's dynamics and optimal solutions. 

It is thus essential to clarify whether the behavior of the EN model observed above changes or persists if 
we relax the simplifying assumption of undistorted retinotopy and allow for retinotopic distortions. By 
analyzing the complete system of equations Eqs. ([3j[4|, we study the EN model exactly as originally 
introduced by Durbin and Mitchison |^. 

We again employ the fact that in the vicinity of a supercritical bifurcation where the non-orientation 
selective state becomes unstable, the entire set of nontrivial fixed points of Eqs. ([sj [4| is determined by the 
third-order terms of the Volterra series representation of the nonlinear operators F^[z,r] and F^[z,r]. The 
model symmetries Eqs. restrict the general form of the leading order terms for any model for the 
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joint optimization of 0PM and RM to 

dtz{^) = L,[z]+Q^[r,z]+7V|[z,z,z] + ... (26) 
a,r(x) = L,[r] + Q^[z,^] + ... . (27) 

Because the uniform retinotopy is linearly stable, retinotopic distortions are exclusively induced by a 
coupling of the RM to the 0PM via the quadratic vector- valued operator Q^[^, z]. These retinotopic 
distortions will in turn alter the dynamics of the 0PM via the quadratic complex- valued operator (5^[r, z]. 
Close to the point of pattern onset (r <C 1), the timescale of 0PM development, r = 1/r, becomes 
arbitrarily large and retinotopic deviations evolve on a much shorter timescale. This separation of 
timescales allows for an adiabatic elimination of the variable r(x), assuming it to always be at the 



equilibrium point of Eq. (27): 



r(x) = -L-MQ'-[^,^]] . (28) 

We remark that as A^^^(/c) < for all finite wave numbers /c > 0, the operator L^[r] is indeed invertible 
when excluding global translations in the set of possible perturbations of the trivial fixed point. Via Eq. 



(28), the coupled dynamics of 0PM and RM is thus reduced to a third-order effective dynamics of the 



0PM: 

dtz{^) ^ L,[z]-Q^[L;'[Q^[z,z]],z]^Ni[z,z,z] 

= L,[z]+N^[z,z,z] + Ni[z,z,z]. (29) 

The nonlinearity N^[z,z,z] accounts for the coupling between 0PM and RM. Its explicit analytical 
calculation for the EN model is rather involved and yields a sum 

12 

N^[z,z,z] = ^aiN^[z,z,z]. 
i=i 

The individual nonlinear operators are nonlinear convolution-type operators and are presented in the 
Methods section together with a detailed description of their derivation. Importantly, it turns out that the 
coefficients are completely independent of the orientation stimulus ensemble. 

The adiabatic elimination of the retinotopic distortions results in an equation for the 0PM (Eq. ([29|)) 



which has the same structure as Eq. (13), the only difference being an additional cubic nonlinearity. Due 
to this similarity, its stationary solutions can be determined by the same methods as presented for the case 
of a fixed retinotopy. Again, via weakly nonlinear analysis we obtain amplitude equations of the form Eq. 
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Figure 8: Angle-dependent interaction functions for the coupling between OPM and RM in the 

EN model, (a, b) griot) and /r(<^) for T]r = 0.005 and a/A = 0.3 (a) and 0.1 (b). (c, d) gr{(^) and fr{(^) 
for r]r = 0.05 and and a/A = 0.3 (c) and 0.1 (d). (e, f) gr{(^) and fr{o^) for r]r = 0.5 and a/A = 0.3 (e) and 
0.1 (f). 



( |15| ). The nonlinear coefficients and fij are determined from the angle-dependent interaction functions 
g{a) and f{a). For the operator N^[z^z,z]^ these functions are given by 

f fl - a^ - 2e-^'^') g2/c>2(cosa-i) ^ g-fc^"^^ 



2a4 [r]r + cr2e-2^c^'(^°««-i)) 



1 



verifying that , N^[z^z^z] is independent of the orientation stimulus ensemble. Besides the interaction 
range a/ A the continuity parameter rjr G [0, oo] for the RM appears as an additional parameter in the 
angle-dependent interaction function. Hence, the phase diagram of the EN model will acquire one 
additional dimension when retinotopic distortions are taken into account. We note, that in the limit 
r]r ^ oo, the functions gr{(^) and fr{(^) tend to zero and as expected one recovers the results presented 
above for fixed uniform retinotopy. The functions gr{ci) and fr{ci) are depicted in Fig. [sjfor various 
interaction ranges a/ A and retinotopic continuity parameters rjr- 
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Figure 9: Coupled n-ECPs as dimension-reducing solutions of the EN model Coupled n-ECP are 
displayed in visual space showing simultaneously the distortion of the retinotopic map and the orientation 
preference map (<j/A = 0.3 (77 = 0.028), 77^ = 7^, circular stimulus ensemble). The distorted grid represents a 
the cortical square array of cells. Each grid intersection is at the receptive field center of the corresponding 
cell. Preferred stimulus orientations are color-coded as in Fig. |2^. As in Fig. [4j n and i enumerate the number 
of nonzero wave vectors and non-equivalent configurations of wave vectors with the same n, respectively. The 
coupled 1-ECP is a pinwheel-free stripe pattern without retinotopic distortion. Only the most anisotropic 
and the most isotropic coupled n-ECPs are shown for each n. Note that for all ECPs, high gradients within 
the orientation mapping coincide with low gradients of the retinotopic mapping and vice versa. Retinotopic 
distortions are displayed on a 5-fold magnified scale for visualization purposes. 
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Coupled Essentially Complex n-Planforms 

In the previous section, we found that by an adiabatic ehmination of the retinotopic distortions in the 



dynamics Eqs. (26, 27) the system of partial integro-differential equations can be reduced to a single 
equation for the 0PM. In this case, the stationary solutions of the 0PM dynamics are again planforms 
composed of a discrete set of Fourier modes 



N 



^(x) = ^A,e^'^^^ (30) 

j 

with |k| = /cc- However, each of these stationary planform 0PM solutions induces a specific pattern of 



retinotopic distortions via Eq. (28). The joint mapping {z(x + r(x),x)} is then an approximate stationary 



solution of Eqs. (26, 27) and will be termed coupled planform solution in the following. In contrast to 



other models for the joint mapping of orientation and visual space (e.g. [3l]33,93j), the coupling between 
the representation of visual space and orientation in the EN model is not induced by model symmetries but 
a mere consequence of the joint optimization of 0PM and RM that requires them to be matched to one 
another. 

For planforms Eq. ( [3Q| ), it is possible to analytically evaluate Eq. ([28| and compute the associated 
retinotopic distortions r(x). After a somewhat lengthy calculation (see Methods), one obtains 

* (^ {AjAk) cos (A^-fcx) + J? {AjAk) sin (A^-fex)) , (31) 



with Ajk = kj — k/e and A£(/c) = —k'^{r]r + e ^ ^ a^). The retinotopic distortions Eq. (31) represent 



superpositions of longitudinal modes (see Fig. Isb). Hence, coupled planform stationary solutions of the EN 



dynamics do not contain any transversal mode components. According to Eq. (31), the pinwheel-free 
coupled 1-ECP state has the functional form {r(x) = 0,z(x) = Aqc'^^^}. This means, that the orientation 
stripe solution does not induce any deviations from the perfect retinotopy as shown previously from 
symmetry. This is not the case for the square pinwheel crystal (sPWC) 

^spwc(x) (X sm{kcXi) + ism{kcX2) , 

the second important solution for undistorted retinotopy. Inserting this ansatz into Eq. (31) and neglecting 
terms of order O ^^e~^^^^^ ^ or higher, we obtain 

rspwc(x) (X / c V c i; 



cr2A5(2fec) V kcSm{2kcX2) 
31 



These retinotopic distortions are a superposition of one longitudinal mode in x-direction and one in 
y-direction, both with doubled wave number ~ 2kc. The doubled wave number implies that the form of 
retinotopic distortions is independent of the topological charge of the pinwheels. Importantly, the gradient 
of the retinotopic mapping R(x) = X + rspwc(x) is reduced at all pinwheel locations. The coupled sPWC 
is therefore in two ways a high coverage mapping as expected. Firstly, the representations of cardinal and 
oblique stimuli (real and imaginary part of z{-k)) are orthogonal to each other. Secondly, the regions of 
highest gradient in the orientation map correspond to low gradient regions in the retinotopic map. 
In Fig. [9| the family of coupled n-ECPs is displayed, showing simultaneously the distortions of the 
retinotopic map and the orientation preference map. Retinotopic distortions are generally weaker for 
anisotropic n-ECPs and stronger for isotropic n-ECPs. However, for all stationary solutions the regions of 
high gradient in the orientation map coincide with low gradient regions (the folds of the grid) in the 
retinotopic map. This is precisely what is generally expected from a dimension-reducing 



mapping 21 62 63,92 . In the following section, we will investigate which of these solutions become 



optimal depending on the two parameters cr/A and r]r that parametrize the model. 
The impact of retinotopic distortions 

According to our analysis, at criticality, the nontrivial stable fixed points of the EN dynamics are 
determined by the continuity parameter r] e (0, 1) for the 0PM or, equivalently, the ratio 
cr/A = ^^/\og(l/r]) and the continuity parameter r]r for the mapping of visual space. We first tested for 
the stability of pinwheel-free orientation stripes (OS) and rhombic pinwheel crystals (rPWCs) solutions of 



Eq. (15), with coupling matrices gij and fij as obtained from the nonlinearities in Eq. (29). The angle 
which minimizes the energy Urpwc (Eq. ([2Q|)) is not affected by the coupling between retinotopic and 
orientation preference map and is thus again a = 7r/4. By numerical evaluation of the criteria for intrinsic 
and extrinsic stability, we found both, OS and sPWC, to be intrinsically and extrinsically stable for all cr/A 
and r]r- 

Next, we tested for the stability of coupled n-ECP solutions for 2 < n < 20. We found all coupled n-ECP 
configurations with n > 2 to be intrinsically unstable for all cr/A and r]r- Evaluating the energy assigned to 
sPWCs and OS identified two different regimes: (i) for shorter interaction range cr/A the sPWC is the 
minimal energy state and (ii) for larger interaction range cr/A the optimum is an OS as indicated by the 
phase diagram in Fig. [TO^. The retinotopic continuity parameter has little influence on the energy of the 
two fixed points. The phase border separating stripes from rhombs runs almost parallel to the ?7^-axis. We 
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Figure 10: Phase diagram of the EN model with variable retinotopy for a circular stimulus 
ensemble 1 57, 64-66, 92] (a) Regions of the ?7r-cr/A-plane in which n-ECPs or rPWCs have minimal 
energy, (b) Pinwheel density time courses for four different simulations of Eqs. ([sj [4| with r = 0.1, 
cr/A = 0.13 (t^ = 0.51), Tjr = Tj (grey traces, individual realizations; red trace, mean value; black trace, 
realization shown in c). (c) OPMs (upper row), their power spectra (middle row), and retinotopic maps 
(lower row) obtained in a simulation of Eqs. ([3j|4|; parameters as in b. (d) Pinwheel density time courses for 
four different simulations of Eqs. (|3j [4| with r = 0.1, a/A = 0.3 (r^ = 0.03), rjr = r] (grey traces, individual 
realizations; red trace, mean value; black trace, realization shown in e). (e) OPMs (upper row), their power 
spectra (middle row), and retinotopic maps (lower row) in a simulation of Eqs. ([3||4|; parameters as in d. 
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numerically confirmed these analytical predictions by extensive simulations of Eq. (|3j [4|) (see Methods for 
details). Fig. \w\ : shows snapshots of a representative simulation with small interaction range (r = 0.1, 
a/A = 0.1{r] = 0.67), r]r = v)- After the initial symmetry breaking phase, the 0PM layout rapidly 
converges towards a crystalline array of pinwheels, the predicted optimum in this parameter regime (Fig. 



lOp). Retinotopic deviations are barely visible. Fig. 10 3 displays pinwheel density time courses for four 
such simulations. Note that in one simulation, the pinwheel density drops to almost zero. In this 
simulation, the OP pattern converges to a stripe-like layout. This is in line with the finding of bistability of 
rhombs and stripes in all parameter regimes. Although the sPWC represents the global minimum in the 
simulated parameter regime, OS are also a stable fixed point and, depending on the initial conditions, may 
arise as the final state of a fraction of the simulations. In the two simulations with pinwheel densities 
around 3.4, patterns at later simulation stages consist of different domains of rhombic pinwheel lattices 
with a < 7r/2. 

Fig. [To]i-e show the corresponding analysis with parameters for larger interaction range r = 0.1, 
cr/A = 0.15 (?7 = 0.41), rjr = V- Here after initial pinwheel creation, pinwheels typically annihilate 
pairwisely and the 0PM converges to an essentially pinwheel-free stripe pattern, the predicted optimal 
solution in this parameter regime (Fig. [lO^). Retinotopic deviations are slightly larger. The behavior of 
the EN model for the joint optimization of RM and 0PM thus appears very similar compared to the fixed 
retinotopy case. Perhaps surprisingly, the coupling of both feature maps has little effect on the stability 
properties of the fixed points and the resulting optimal solutions. 



As in the previous case, the structure of the phase diagram in Fig. 10 i appears somewhat counterintuitive. 
A high coverage and pinwheel-rich solution is the optimum in a regime with large 0PM continuity 
parameter where discontinuities in the 0PM such as pinwheels should be strongly penalized. A 
pinwheel-free solution with low coverage and high continuity is the optimum in a regime with small 
continuity parameter. As explained above, a large 0PM continuity parameter at pattern formation 



threshold implies a small interaction range a/ A (see Eq. (24)). In such a regime, the gain in coverage by 
representing many orientation stimuli in a small area spanning the typical interaction range, e.g. with a 
pinwheel, is very high. Apparently this gain in coverage by a regular positioning of pinwheels outweighs 
the accompanied loss in continuity for very large 0PM continuity parameters. This counterintuitive 
interplay between coverage and continuity thus seems to be almost independent of the choice of retinotopic 
continuity parameters. 

The circular orientation stimulus ensemble contains only stimuli with a fixed and finite "orientation 
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energy" or elongation \sz\. This raises the question of whether the simple nature of the circular stimulus 
ensemble might restrain the dynamics of the EN model. The EN dynamics are expected to depend on the 
characteristics of the activity patterns evoked by the stimuli and these will be more diverse and complex 
with ensembles containing a greater diversity of stimuli. Therefore, we repeated the above analysis of the 
EN model for a richer stimulus ensemble where orientation stimuli are uniformly distributed on the disk 
{sz, \sz\ < 2}, a choice adopted by a subset of previous studies, e.g. p^[25|[82|[94| . In particular, this 
ensemble contains unoriented stimuli with l^^l =0. Intuitively, the presence of these unoriented stimuli 
might be expected to change the role of pinwheels in the optimal 0PM layout. Pinwheels' population 
activity is untuned for orientation. Pinwheel centers may therefore acquire a key role for the representation 
of unoriented stimuli. Nevertheless, we found the behavior of the EN model when considered with this 
richer stimulus ensemble to be virtually indistinguishable from the circular stimulus ensemble. Details of 
the derivations, phase diagrams and numerically obtained solutions are given in Appendix I. 



Are there stimulus ensembles for which realistic, aperiodic maps are optimal? 

So far, we have presented a comprehensive analysis of optimal dimension-reducing mappings of the EN 
model for two widely used orientation stimulus distributions (previous sections and Appendix I). In both 
cases optima were either regular crystalline pinwheel lattices or pinwheel-free orientation stripes. These 
results might indicate that the EN model for the joint optimization of 0PM and RM is per se incapable of 
reproducing the structure of orientation preference maps as found in the visual cortex. Drawing such a 
conclusion is suggested in view of the apparent insensitivity of the model's optima to the choice of stimulus 
ensemble. The two stimulus ensembles considered so far however do not exhaust the infinite space of 
stimulus distributions that are admissible in principle. From the viewpoint of "biological plausibility" it is 
certainly not obvious that one should strive to examine stimulus distributions very different from these, as 
long as the guiding hypothesis is that the functional architecture of the primary visual cortex optimizes the 
joint representation of the classical elementary stimulus features. If, however, stimulus ensembles were to 
exist, for which optimal EN mappings truly resemble the biological architecture, their characteristics may 
reveal essential ingredients of alternative optimization models for visual cortical architecture. 
Adopting this perspective raises the technical question of whether an unbiased search of the infinite space 
of stimulus ensembles only constrained by the model's symmetries Eqs. is possible. To answer this 



question, we examined whether the amplitude equations Eq. (15) can be obtained for an arbitrary 
orientation stimulus distribution. Fortunately, we found that the coefficients of the amplitude equations are 
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completely determined by the finite set of moments of order less than 5 of the distributions. The approach 
developed so far can thus be used to comprehensively examine the nature of EN optima resulting for any 
stimulus distribution with finite fourth-order moment. While such a study does not completely exhaust the 
infinite space of all eligible distributions, it appears to only exclude ensembles with really exceptional 
properties. These are probability distributions with diverging fourth moment, i.e. ensembles that exhibit a 
heavy tail of essentially "infinite" orientation energy stimuli. 

Since the coupling between 0PM and RM did not have a large impact in the case of the two classical 
stimulus ensembles, we start the search through the space of orientation stimulus ensembles by considering 
the EN model with fixed retinotopy r(x) = 0. The coefficients for the nonlinear operators Nf[z^ z] in 



Eq. (14) for arbitrary stimulus ensembles are given by 



_ {\sA') {\s.\') , 1 _ {\sA') {\s.\') {\sA') , {\s.\') 

"1 ~ 16cr6 2cr4 + 2cr2 ^2 — g^^e 327rcr8 "3 — 647rcr8 ^ 167rcr6 

"7 - 487r2aio 247r2a8 - ^Qt,'2^io "9 - 2567r3ai2 

_ {\s.\') {\s.\^) _ (k.l^) 

"10 ~ 487r2crio 247r2cr8 "H ~ 967r2crio * 

The corresponding angle-dependent interaction functions are given by (see Methods) 



(32) 



2(j4 



|2\ 

fi_Z 1^1 _ e-2'=c<^' (cosh(2fc2o-2cosa) +2cosh(fc^(T2cosa)) +2e-'=''" j 
+ ^ (e-2'==-' cosh(2fc2cT2 cos a) - l) + i^e'^'^'"' sinh* {ll2kla^ cos a) . (33) 

Again, without loss of generality, we set (l^^p) = 2. At criticality, both functions are parametrized by the 



continuity parameter 77 G (0, 1) for the 0PM or, equivalently, the interaction range a/A = ^/\og{)Jr]) 
and the fourth moment (l^^l^) of the orientation stimulus ensemble. The fourth moment, is a measure of 
the peakedness of a stimulus distribution. High values generally indicate a strongly peaked distribution 
with a large fraction of non-oriented stimuli {\sz\^ ~ 0), together with a large fraction of high orientation 
energy stimuli {\sz\^ large). 

The dependence of g{a) on the fourth moment of the orientation stimulus distribution and f{a) suggests 
that different stimulus distributions may indeed lead to different optimal dimension-reducing mappings. 
The circular stimulus ensemble possesses the minimal possible fourth moment, with 
(I'^^l^) = ((k^P))^ = 4. The fourth moment of the uniform stimulus ensemble is (l^^l^) = 16/3. The 
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Figure 11: Angle-dependent interaction functions for the EN model with fixed retinotopy for dif- 
ferent fourth- moment values of the orientation stimulus distribution and effective interaction- 
widths, (a, b) g{a) and f{a) for S4 = 8 and a/A = 0.1 (a) and a/A = 0.3 (b). (c, d) g{a) and f{a) for 
S4 = 20 and a/A = 0.1 (c) and a/ A = 0.3 (d). (e, f) g{a) and f{a) for 54 = 100 and a/A = 0.1 (e) and 
cr/A = 0.3 (f). 
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angle-dependent interaction functions for both ensembles (Eq. (25), Fig. 22 in Appendix I) are recovered, 



when inserting these values into Eq. (33). 



To simplify notation in the following, we define 

S4 = (kH-(|s.P)' = (|s,0-4 
as the parameter characterizing an orientation stimulus distribution. This parameter ranges from zero for 



the circular stimulus ensemble to infinity for ensembles with diverging fourth moments. Fig. 11 displays 
the angle-dependent interaction functions for different values of a/A and 54. In all parameter regimes, g{a) 
and f{a) are larger than zero. The amplitude dynamics are therefore guaranteed to converge to a stable 
stationary fixed point and the bifurcation from the nonselective fixed point in the EN model is predicted to 
be supercritical in general. 

By evaluating the energy assigned to the rPWC and n-ECPs, we investigated the structure of the 
two-dimensional phase space of the EN model with an arbitrary orientation stimulus distribution. Firstly, 
it is not difficult to show that the angle a which minimizes the energy UrPWC (Fq. ( |2Q[ )) of a rhombic 
pinwheel crystal is a = 7r/4 for all a/ A and 54. Hence, a square lattice of pinwheels is in all parameter 



regimes energetically favored over any other rhombic lattice configuration of pinwheels. Fig. 12 displays 
the phase diagram of the EN model with an arbitrary orientation stimulus distribution. For orientation 
stimulus distributions with small fourth moments, optimal mappings consist of either parallel pinwheel-free 
stripes or quadratic pinwheel crystals. These distributions include the circular and the uniform stimulus 
ensembles with S4 = and 54 = 4/3. Above a certain value of the fourth moment around 54 = 6, n-ECPs 
with n > 2 become optimal mappings. For a short interaction range cr/A, hexagonal pinwheel crystals 
dominate the phase diagram in a large region of parameter space. With increasing interaction range, we 
observe a sequence of phase transitions by which higher n-ECPs become optimal. For n > 3, these optima 
are spatially aperiodic. In all parameter regimes, we found that the n-ECP with the most anisotropic mode 
configuration (Fig. [4]3, left column) is the energetically favored state for n > 3. Pinwheel densities of these 
planforms are indicated in Fig. [12] and are typically smaller than 2.0. We note that this is well below 
experimentally observed pinwheel density values 3^. 

Optimal mappings of orientation preference for finite fourth moment in the EN model are thus either 
orientation stripes, periodic arrays of pinwheels (hexagonal, square), or aperiodic pinwheel arrangements 
with low pinwheel density. 

We numerically tested these analytical predictions by simulations of Eq. (|3| (r(x) = 0) with two additional 
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0.1 0.2 0.3 0.4 0.5 0.6 q/A 

Figure 12: Stripe-like, crystalline, and quasi-crystalline cortical representations as optimal solu- 
tions to the mapping of orientation preference with fixed uniform retinotopy in the EN model. 

The graph shows the regions of the 54-cr/A-plane in which n-ECPs or sPWCs have minimal energy. For 
n > 3, pinwheel densities of the energeticahy favored n-ECP configuration are indicated. 
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stimulus ensembles with 54 = 6 and 54 = 8 (see Methods). Fig. 13 i shows snapshots of a simulation with 



(r = 0.1, cr/A = 0.2 (7^ = 0.2)) and 54 = 6. After the initial phase of pattern emergence, the 0PM layout 



converges towards an arrangement of fractured stripes which resembles the 2-ECP state (Fig. 13 i, most 
right), the optimum predicted in this regime. In the Fourier spectra, two distinct peaks of the active modes 
are clearly visible in the final stages of the simulation (Fig. p!3^ , lower row). The 2-ECP state is exotic in 
the sense that it is the only n-ECP containing line defects and thus the pinwheel density is not a 
well-defined quantity. This explains the pronounced numerical variability in the measured pinwheel 
densities in simulations during the convergence towards a 2-ECP state (Fig. \i3\p). 

Fig. \i3\ : shows snapshots of a simulation with (r = 0.1, cr/A = 0.2(7^ = 0.2)) and 54 = 8 (Gaussian stimulus 
ensemble). After the initial phase of pattern emergence, the 0PM layout converges towards a regular 



hexagonal arrangement of pinwheels which resembles the anisotropic 3-ECP (Fig. 13 3, far right), the 



optimum predicted in this regime. In the Fourier spectra, three distinct peaks forming an angle of 60 



degrees are clearly visible in the later stages of the simulation (Fig. 13 3, lower row). Pinwheel densities in 



the simulations consistently approach the theoretically predicted value of 2cos(7r/6) 1.73 (Fig. [Tsji). 
Permutation symmetric limit 

In the previous section, we uncovered a parameter regime for the EN model in which optimal solutions are 
spatially aperiodic. This can be viewed as a first step towards realistic optimal solutions. In the identified 
regime, however, among the family of n-ECPs only those with pinwheel densities well below experimentally 



observed values 38 are energetically favored (see Fig. 12). In this respect, the repertoire of aperiodic 



optima of the EN model differs from previously considered abstract variational models for 0PM 



development [35y36[|38[|39| . In these models, an energetic degeneracy of aperiodic states with low and high 
pinwheel densities has been found which leads to a pinwheel statistics of the repertoire of optimal solutions 
that quantitatively reproduces experimental observations |38j|95]. What is the reason for this difference 
between the two models? In 35 , the energetic degeneracy of aperiodic states with low and high pinwheel 
densities was derived from a so-called permutation symmetry 

Ni[u,v,w]=N^[w,u,v], (34) 

of the cubic nonlinearities of the model. It can be easily seen, that the cubic nonlinearities obtained in the 
third order expansion of the EN model do not exhibit this permutation symmetry (see Methods). As 
shown by Reichl 96 , the absence of permutation symmetry can lead to a selection of a subrange of 



40 



a t=1x t=10T t = 50T t = 600x 




time [x] 



Figure 13: Approaching crystalline n-ECP optima in the EN model with fixed retinotopy. (a) 

OPMs (upper row) and their power spectra (lower row) in a simulation of Eq. ([3| with r(x) = 0, r = 0.1, 
cr/A = 0.2 and = 6. The predicted optimum is the 2-ECP (black frame), (b) Pinwheel density time courses 
for four different simulations (parameters as in a; gray traces, individual realizations; black trace, simulation 
in a; red trace, mean value) (c) OPMs (upper row) and their power spectra (lower row) in a simulation of 
Eq. ([3| with r(x) = 0, r = 0.1, a/ A = 0.3 and 54 = 8. The predicted optimum is the anisotropic 3-ECP 
(black frame), (d) Pinwheel density time courses for four different simulations (parameters as in c; gray 
traces, individual realizations; black trace, simulation in c; red trace, mean value). 
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Figure 14: Quantifying permutation symmetry breaking in the EN model, (a-d) g{a) (red traces) 
and the "permutation symmetrized" function gpm{(^) = l/2{g{a) -\- g{a -\- n)) (blue traces, see Eq. (35)) for 
cr/A = 0.1 and 0.3 and 54 = 6 and 100. (e) Permutation symmetry parameter d (Eq. (36)) in the EN model 
with fixed retinotopy. Permutation symmetry breaking is largest for a/ A ^ 0.25 and small 54. In the limit 
54 ^ 00, permutation symmetry is restored. 



pinwheel densities in the repertoire of optima of 0PM models. Depending on the degree of permutation 
symmetry breaking, the family of optima of such models, albeit encompassing aperiodic 0PM layouts, may 
consist of layouts with either unrealistically low or high pinwheel densities. Furthermore, for very strong 
permutation symmetry breaking, stationary solutions from solution classes other than the n-ECPs and 
rPWCs with low or high pinwheel densities may become optima of models for 0PM development. In order 
to determine a regime in which the EN model optima quantitatively resemble experimentally observed 
0PM layouts, it is therefore important to quantify the degree of permutation symmetry breaking in the 
EN model and to examine whether permutation symmetric limits exist. As shown in the Methods section. 



any cubic nonlinearity N^[z^z^z] that obeys the permutation symmetry Eq. (34) has a corresponding 
angle-dependent interaction function g{a) which is 7r-periodic. Therefore, we examine the degree of 
permutation symmetry breaking in the EN model by comparing the angle-dependent interaction function 



g{a) of its third order expansion (see Eq. (33) and Fig. 11) to the 7r-periodic function 

9pm{oi) = 1/2 {g{a) + ^(a + tt)). This "permutation- symmetrized" part of the angle-dependent interaction 
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function of the EN model for general orientation stimulus ensembles reads 

g^M = ^^^^e-2^'^%inli^(l/2i^,V2cosa) 

-i\^e-^^>" ((cosh {2eycosa) - 2cosh (i^.V^cosa)) - 2e^>" - e^^>"^ 

+ ^ (l + e"^^'^' cosh(2/c^cr^ cosa)^ . (35) 



A comparison between gpm{o^) and g{a) is depicted in Fig. 14i-d. It shows that essentially insensitive to 
the interaction range cr/A, at large values of the fourth moment original and permutation symmetrized 
angle-dependent interaction functions converge. We quantified the degree of permutation symmetry 
breaking with the parameter 

d= ll%^pl^ gn(^(Q)-^(^)). (36) 

This parameter is zero in the case of a permutation symmetric cubic nonlinearity. In the case of a 
g- function completely antisymmetric around a = 7r/2, the parameter is either plus or minus one, depending 
on whether the maximum of gpm is at zero or tt. If d is smaller than zero, low pinwheel densities are 
expected to be energetically favored and vice versa. The values of d in parameter space is depicted in Fig. 



143. It is smaller than zero in the entire phase space, implying a tendency for low pinwheel density optimal 



states, in agreement with the phase diagram in Fig. 12 Permutation symmetry breaking is largest for cr/A 
around 0.25 and small fourth moment values of the orientation stimulus distribution. It decays to zero for 



large fourth moments proportionally to I/54 as can be seen by inserting Eqs. (33) and (35) into Eq. (36) . 
In the infinite fourth moment limit 54 ^ 00, the cubic nonlinearities of the third order expansion of the EN 
model become permutation symmetric. 

In this case, the EN model is parametrized by only one parameter, the effective intracortical interaction 



range a/ A and we obtain a rather simple phase diagram (Fig. 15). Optimal solutions are n-ECPs for 
increasing a/ A and we observe a sequence of phase transitions towards a higher number of active modes 
and therefore more complex spatially aperiodic 0PM layouts. Importantly, for a subregion in the phase 
diagram with given number of active modes, all possible n-ECP mode configurations are energetically 
degenerate. It is precisely this degeneracy that has been previously shown to result in a pinwheel statistics 
of the repertoire of aperiodic optima which quantitatively agrees with experimental observations 38 . 
Therefore, our unbiased search in fact identified a regime, namely a very large effective interaction range 
and infinite fourth moment of the orientation stimulus ensemble, in which the EN model formally predicts 
which quantitatively reproduce the experimentally observed VI architecture. 
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Figure 15: Phase diagram of the EN model with fixed retinotopy in the permutation symmetric 
hmit 54 ^ oo. The graphs show the regions on the a/ A- axis (lower axis) and the corresponding 77- axis 
(upper axis), where n-ECPs or sPWCs have minimal energy. High n-ECPs (n > 10) exhibit universal 
pinwheel statistics. Note however the extremely small 77- values for large <j/A. 

Unexpectedly, however, this regime coincides with the limit of applicability of our approach. Permutation 
symmetry is exactly obtained by approaching stimulus distribution with diverging fourth moment for 
which the amplitude equations may become meaningless. We would generally expect that the EN for very 
large but finite fourth moment can closely resemble a permutation symmetric model. However, to 
consolidate the relevance of this regime, it appears crucial to establish the robustness of the limiting 
behavior to inclusion of retinotopic distortions. 



Optimal solutions of the EN model with variable retinotopy and arbitrary orientation stimulus 
ensembles 

In the EN model for the joint mapping of visual space and orientation preferences, the angle-dependent 
interaction functions depend on four parameters: 7^, a, the fourth moment (l^^l^) of the stimulus ensemble 
and r]r' By setting a = (J*{r])^ we are left with three free parameters at criticality. Therefore, a 
three-dimensional phase diagram now completely describes pattern selection in the EN model. For better 
visualization, in Fig. [16] we show representative cross sections through this three-dimensional parameter 
space for fixed rfr. Firstly, we note the strong similarity between the phase diagram for fixed retinotopy 



(Fig. 12) and the cross sections through the phase diagrams for the joint mappings shown in Fig. 16 This 
expresses the fact that retinotopic mapping and 0PM are only weakly coupled or mathematically, 
gr{oi) <C g{a) in all parameter regimes (see Appendix HI). Again, for distributions with small fourth 
moment, optimal mappings consist of either pinwheel-free orientation stripes or sPWCs. Above a certain 
fourth moment value around S4 = 6, higher coupled n-ECPs are optimal. For small interaction range a/ A, 
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Figure 16: Stripe-like, crystalline, and quasi-crystalline cortical representations as optimal solu- 
tions to the joint mapping problem of visual space and orientation preference in the EN. (a-d) 

Phase diagrams for the joint mapping of visual space and orientation preference in the EN near criticahty 
for 77^ = (a), r]r = 0.01 (b), r]r = 0.1 (c), and rjr = 10 (d). The graphs show the regions of the 54-cr/A-plane 
in which coupled n-ECPs or sPWCs have minimal energy. For n > 3, pinwheel densities of the energetically 
favored n-ECP configuration are indicated. Note the strong similarity between the phase diagrams and the 
phase diagrams in the fixed retinotopy case (Fig. 12). 
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hexagonal pinwheel crystals (coupled 3-ECPs) represent optimal mappings in a large fraction of parameter 
space. With increasing a/A, we observe a sequence of phase transitions by which higher n-ECPs become 
optimal. Anisotropic planforms at the lower end of the spectrum of pinwheel densities are always 
energetically favored over high pinwheel density layouts. The only difference between the cross-sections is 
that the region covered by sPWCs increases for decreasing rjr- The phase diagram for large r]r = 10 is 



virtually indistinguishable from the phase diagram in Fig. 16 

Optimal mappings of orientation preference are thus either orientation stripes, periodic arrays of pinwheels 
(hexagonal, quadratic) or quasiperiodic pinwheel arrays with low pinwheel density. Retinotopic distortions 
lead to lower gradients of the retinotopic mapping at high gradient regions of the 0PM. This is in line with 



some of the experimental evidence [55y97| but contradicts others [98 . 

Most importantly, we note that the results on permutation symmetry breaking in the fixed retinotopy case 
are not altered by allowing for retinotopic distortions. Since gr{(^) does not depend on the fourth moment 
of the orientation stimulus distribution, non-permutation symmetric terms decay as l/s^ for large s^. 



Hence, in the limit 54 ^ oo, permutation symmetry is restored and we recover the phase diagram Fig. 15 
also for the EN model with variable retinotopy independent of rjr- As the energy contribution of 
retinotopic deviations r(x) becomes negligible in the infinite fourth moment limit, the optima are then 
simply the corresponding coupled n-ECPs and these states are energetically degenerate for fixed n. For 
very large effective interaction range and infinite fourth moment of the orientation stimulus ensemble, the 
EN model with variable retinotopy is able to quantitatively reproduce the experimentally observed 
pinwheel statistics in OPMs. It furthermore predicts reduced gradients of the visual space mapping at high 
gradient regions of the 0PM. 



Finite stimulus samples and discrete stimulus ensembles 

Our reexamination of the EN model for the joint optimization of position and orientation selectivity has 
been so far carried out without addressing the apparently fundamental discrepancy between our results 
and the large majority of previous reports. Since the seminal publication of Durbin and Mitchison in 
1990 [25, numerous studies have used the EN model to simulate the development of visual cortical maps or 
to examine the structure of optimal mappings by numerical simulation 58,62-65,80,99j[l00 . These studies 



have either used the circular or the uniform orientation stimulus ensemble for which, to the best of our 
knowledge, the only two nontrivial stationary solutions are square pinwheel crystals or orientation stripes. 
Furthermore, we found that the gradient descent dynamics seems to readily converge to the respective 
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a 216 stimuli |j 10^ stimuli 




Figure 17: Development of OPM and retinotopic distortions in EN simulations with fixed 
stimulus sets of different sizes, (a) OPMs (left), their power spectra (middle) and retinotopic maps 
(right) for t = lOr (upper row) and t = 300r (lower row) obtained in simulations with fixed stimulus set 
[t] = 0.028, cr/A = 0.3, s4 = 4/3, 216 stimuh). (b) 10^ stimuh (all other parameters as in a), (c) 10^ stimuh 
(all other parameters as in a). Large stripe-like OP domains are generated via pairwise pinwheel annihilation 
for large simulation times. Retinotopic distortions are fairly weak.(d) Pinwheel density time course for EN 
simulations with fixed stimulus sets of different sizes, including the simulations from a-c (red, green, blue 
traces 216, 10^, 10^ stimuli) (all other parameters as in a). Dashed lines represent individual simulations, 
solid lines an average over four simulations. Note, that the pinwheel density rapidly decays below 2.0 in 
both cases, and in particular for 10^ stimuli, the OPM pattern acquires large stripe-like regions. 
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minima of the EN free energy. This indicates that other local minima and more complex intrinsically 
aperiodic states are not dominant in this model. In fact, we found that all aperiodic stationary solutions 
we could perturbative calculate analytically are unstable and thus represent hyperbolic saddle points and 
not local minima (see also Fig. [21^,b). As these stable solutions barely resemble experimentally observed 
orientation preference maps, it is not obvious how the EN model in all of these studies could appear as a 
model well-suited to describe the complex layout of real cortical orientation maps. Prior studies however 
often used computational methods different from our fixed parameter steepest descent simulations. 
Two alternative approaches have been used predominantly to study dimension reducing mappings for 
cortical representations. These methods have been applied to both the EN model and the other widely 
used dimension reduction model, the Self-Organizing Feature Map (SOFM), originally introduced by 
Kohonen 59 . The simplest way to compute mappings from a high dimensional feature space onto the 2D 
model cortex is by iterating the following procedure for a large number of randomly chosen stimuli 
(e.g. [56l[57l[66H68l [TM|[TQ2] ) : (i) Stimuh are chosen one at a time randomly from the complete feature 
space, (ii) The activation function for a particular stimulus is computed. In the case of the EN model, this 
activation function can acquire a rather complex form with multiple peaks (see Discussion). In the case of 
an SOFM, this activation function is a 2D-Gaussian. (iii) The preferred features of the cortical grid points 
are updated according to a discretized version of Eqs. (|3j [4| or the corresponding equations for the SOFM 
model. Typically, this procedure is repeated on the order of 10^ times. The resulting layout is then 
assumed to at least approximately solve the dimension reduction problem. In many studies, small stimulus 
sets have been chosen presumably for computational efficiency and not assuming specifically that the 



cortex is optimized for a discrete finite set of stimuli. In 21 for instance, a set of 216 stimuli was used, 
that was likely already at the limit of computing power available at this time. 

In a more refined approach, the EN model as well as Kohonen's SOFM model have been trained with a 
finite set of stimuli (typically with on the order of 10^-10^) and the final layout of the model map has been 
obtained by deterministic annealing 103 , i.e. by gradually reducing the numerical value of a in a 



numerical minimization procedure for the energy functional T at each value of a (see e.g. 21 ,64 , 65[|8Q| and 
Methods). In such simulations, often non-periodic boundary conditions were used. One might suspect in 
particular the second approach to converge to 0PM layouts deviating from our results. It is conceivable in 
principle, that deterministic annealing might track stationary solutions across parameter space that are 
systematically missed by both, our continuum limit analytical calculations as well as our descent numerical 
simulations. 
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To assess the potential biases of the different approaches, we implemented (i) finite stimulus sampling in 
our gradient descent simulations and (ii) studied the results of deterministic annealing simulations varying 
both the size of the stimulus set as well as the type of boundary conditions applied. 
We simulated the dynamics Eqs. ([3)[4]) with finite sets of stimuli of different sizes (see Methods), drawn 



from the circular stimulus ensemble. As e.g. in 21,25 , was set to a small value {r] = 0.025) such that 
the optimal configuration for the joint mapping of visual space and orientation preference is the coupled 



1-ECP (see Fig. 10), i.e. a pattern of parallel orientation stripes without any retinotopic distortion (see 
Fig. [9|. Fig. 17 displays representative simulations for stimulus sets of size N = 216 (as used in [21 ) (a), N 
= 10^ (b), N = 10^ (c) stimuli. Simulation time t is measured in units of the intrinsic time scale r (see 
Methods). For N = 216 stimuli, RM and 0PM quickly reach an apparently stationary configuration with a 
large number of pinwheels at around t = 20r. Power is distributed roughly isotropically around the origin 
of Fourier space (k = 0). The stable 0PM lacks a typical length scale and, expressing the same fact, the 
power spectrum lacks the characteristic ring of enhanced Fourier amplitude. Retinotopic distortions are 
fairly pronounced. Both obtained maps resemble the configurations reported in 21 . 

For N = 10^ stimuli, we find that OPMs exhibit a characteristic scale (see dark shaded ring in the power 
spectrum) and a dynamic rearrangement of the maps persists at least until t = 200r. Stripe-like OP 
domains are rapidly generated via pairwise pinwheel annihilation for t > lOr. Retinotopic distortions are 
fairly weak. For N = 10^ stimuli, again OPMs exhibit a characteristic scale (see dark shaded ring in the 
power spectrum) and the map dynamics persists beyond t = 200r. A larger fraction of the pinwheels 
annihilate pairwisely compared to N = 10^ stimuli, leading progressively to a pattern with large stripe-like 
domains. Retinotopic distortions are fairly weak. For both cases with massive stimulus sampling {N = 10^, 
N = 10^), the pinwheel density rapidly drops below the range observed in tree shrews, galagos and ferrets 
and than further decreases during subsequent map rearrangement. 

In summary, the more stimuli are chosen for the optimization procedure, the less pinwheels are preserved 
in the pattern of orientation preference and the more the resulting map resembles the analytically obtained 
optimal solution. 

Deterministic annealing approaches which change parameters of the energy functional during the 
computational minimization process differ more fundamentally from our gradient descent simulations than 
the iterative schemes used with fixed parameters. Studies using deterministic annealing in addition 



frequently used non-periodic boundary conditions (e.g. 64,65 80]). To study all potential sources of 
deviating results, we implemented deterministic annealing for the EN energy function (see Methods, Eq. 
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46) for periodic boundaries, non-periodic boundary conditions as well as random and grid-like finite 



stimulus ensembles (see Methods). We closely follow the refined methods used in 64 65,80 and performed 
deterministic annealing simulations for the EN model with retinotopic distortions and stimuli drawn from 
the circular stimulus ensemble. 

Figs. 18 i and 19 i display representative simulations for random stimulus sets of size N = 10^, N = 10^ 



and N = 10^ for periodic boundary conditions (Figs. 18 i) and non-periodic boundary conditions (Figs. 
19^). Furthermore depicted are the pinwheel densities of stationary solutions as well as their energies, 
relative to the energy of a pinwheel-free stripe solution (see Methods) for different annealing rates ^ (Figs. 



18 D-d, Figs. |l9p-d). Figs. 183-g and 193-g additionally show the statistics of nearest neighbor (NN) 



pinwheel distances as well as the standard deviation (SD) of the pinwheel densities for randomly selected 
subregions in the 0PM as introduced in 38 , averaged over four simulations with N = 10^. To facilitate 



comparison, we superimposed fits to the experimentally observed statistics ^ for orientation maps in tree 
shrews, ferrets and galagos. 

When annealing with periodic boundary conditions, the maps found with deterministic annealing 
essentially resemble our gradient descent dynamics simulations. The larger the set of stimuli, the more 
stripe-like are the orientation preference maps obtained (Fig. fl8k,b). Furthermore, the more carefully we 



annealed, the lower the pinwheel density of the obtained layouts (Fig. 18 3). For N = 10^, the pinwheel 
density averaged over four simulations with annealing rate 0.999 was p = 2.04 As expected, the energy of 



the final layouts decreased with slower annealing rates (Fig. 181). However, when starting from random 
initial conditions, the energy of the final layouts found was always higher compared to the energy of a 
pinwheel-free stripe solution (see Methods for details), which is the predicted optimum for the circular 
stimulus ensemble. NN-pinwheel distance histograms are concentrated around half the typical column 
spacing and in particular pinwheel pairs with short distances are lacking completely (Fig. [l8^,f). 
For non-periodic boundary conditions and random stimuli, we found that retinotopic distortions are more 
pronounced than for periodic boundary conditions. They however decreased with increasing number of 
stimuli. For large the stimulus numbers, we observed stripe-like orientation preference domains which are 
interspersed with lattice-like pinwheel arrangements (see Fig. [l9]3), lower row, upper left corner of the 
0PM). For N = 10^, the pinwheel density averaged over four simulations with annealing rate 0.999 was 
p = 2.71. 

Similarly to the results for periodic boundary conditions, short distance pinwheel pairs occur less 
frequently than in the experimentally observed maps, indicating an increased regularity in the pinwheel 
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distances compared to real OPMs (Fig. 193,f). This regularity is further indicated by a smaller exponent of 
the SD compared to the Poisson process (Fig. [l9^). The perfect stripe-like solution is not the optimum for 
non-periodic boundaries. The energy of the map layouts found with very slow annealing rates is slightly 
lower than the energy of the pinwheel-free 0PM layout (Figure [l9]i). We note that the layout of the 0PM 
at the boundaries does not differ substantially from the layout inside the simulated domain, suggesting 
that boundary effects affect the entire simulated domain for the relatively small region treated. 
Finally, we performed simulations with grid-like stimulus patterns as e.g. used in f64,65 . These 
simulations displayed a strong tendency towards rhombic pinwheel arrangements, i.e. the second stable 
stationary solution found for the circular stimulus ensemble. We refer to Appendix II for further details. 
In summary, our results for the discrete EN model with deterministic annealing largely agree with the 
analytical results. Irrespective of the numerical methodology, the emerging map structure for large 
numbers of stimuli is confined to the states predicted by our analytical treatment of the continuum 
formulation of the EN. This behavior is expected because the energies underlying the deterministic 
annealing and the steepest descent simulations are mathematically equivalent (see Methods). In any kind 
of deterministic annealing simulation we tested, resulting patterns were patchworks of the two fundamental 
stable solutions identified by the analytical treatment: pinwheel free stripes and square lattices of 
pinwheels. Such patchworks are spatially more complicated than perfect stripes or crystals. Nevertheless, 
they qualitatively differ in numerous respects from the experimentally observed spatial arrangements (see 



Figs. 18, 19 and 27 in Appendix II). How the fundamental stable solutions are stitched together somewhat 
differs between the different kinds of simulations. For instance, using a grid-like stimulus ensemble with 
non-periodic boundary conditions apparently energetically favors the rhombic pinwheel crystal compared 



to the pinwheel-free stripe regions (see Fig. 26 in Appendix II). In summary, while some of the patterns 
obtained by deterministic annealing might be called "good-looking" maps, all of them substantially deviate 
from the characteristics of experimentally observed pinwheel arrangements. 

We conclude, that the differences between our results and those of previous studies are most likely due to 
the small finite stimulus samples used largely for reasons of computational tractability. Deterministic 
annealing using stimulus samples that fill the feature space converges to the same types of patterns found 
by perturbation theory. We further conclude, that our methods do not systematically miss biologically 
relevant local minima of the classical EN energy function. 
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Periodic boundary conditions 




Figure 18: The EN model with periodic boundary conditions, solved with deterministic an- 
nealing (a) OPMs (left) and retinotopic maps (right) for = 10^ (upper row), N = 10^ (middle row) 
and N = 10^ (lower row) random stimuli and periodic boundary conditions (annealing rate x = 0.999, see 
Methods). P is the continuity parameter in the conventional definition of the EN model (see Methods, Eq. 



46) and is scaled, such that a comparable number of columns is emerging in the simulations for each size 



of the stimulus set. (b) Pinwheel densities of EN solutions for different numbers of stimuli, x = 0.999. (c) 
Pinwheel densities of EN solutions for 10^ stimuli and different annealing rates, (d) Energies of solutions 
for 10^ stimuli, relative to the energy of a pinwheel-free stripe solution (see Methods) for different annealing 
rates, (b-d) Crosses mark individual simulations, red line indicates average values, (e-f) Statistics of nearest 
neighbor pinwheel distances for pinwheels of (e) arbitrary and (f) opposite and equal charge for 10^ random 
stimuli and periodic boundary conditions, averaged over four simulations (red curves). Black curves repre- 
sent fits to the experimental data from [38 . (g) Standard deviations (SD) of pinwheel densities estimated 
from randomly selected regions in the 0PM. Black dashed curve indicates SD for a two-dimensional Poisson 
process of equal density. 
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Non-periodic boundary conditions 



N = 103 stimuli, (3 = 2 
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Figure 19: The EN model with non-periodic boundary conditions, solved with deterministic 



annealing (a-g) As Fig. 18, but for non-periodic boundary conditions. 
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Discussion 
Summary 

In this study, we examined the solutions of what is perhaps the most prominent optimization model for the 
spatial layout of orientation and retinotopic maps in the primary visual cortex, the Elastic Network (EN) 
model. We presented an analytical framework that enables us to derive closed-form expressions for 
hyperbolic fixed points, local and global minima, and to analyze their stability properties for arbitrary 
optimization models for the spatial layout of OPMs and RMs. Using this framework, we systematically 
reexamined previously used instantiations of this model, dissecting the impact of stimulus ensembles and of 
interactions between the two maps on optimal map layouts. To our surprise, the analysis yielded virtually 
identical results for all of these model instantiations that substantially deviate from previous numerical 
reports. Pinwheel-free orientation stripes and crystalline square lattices of pinwheels are the only optimal 
dimension-reducing 0PM layouts of the EN model. Both states are generally stable but exchange their 
roles as optima and local minima at a phase border. Numerical simulations of the EN gradient descent 
dynamics as well as simulations utilizing deterministic annealing confirmed our analytical results. For both 
processes, the initially spatially irregular layouts rapidly decayed into a patchwork of stripe- like or 
crystal-like local regions that then became globally more coherent on longer timescales. Pinwheel-free 
solutions were approached after an initial phase of pattern emergence by pairwise pinwheel annihilation. 
Crystalline configurations were reached by the generation of additional pinwheels and pinwheel annihilation 
together with a coordinated rearrangement towards a square lattice. These results indicate that layouts 
which represent an optimal compromise of coverage and continuity for retinotopy and orientation do not 
per se reproduce the spatially aperiodic and complex structure of orientation maps in the visual cortex. 
To clarify whether the EN model is in principle capable of reproducing the biological observations, we 
performed an unbiased comprehensive inspection of EN optima for arbitrary stimulus distributions 
possessing finite fourth moments. This analysis identified two key parameters determining pattern 
selection: (i) the effective intracortical interaction range and (ii) the fourth moment of the orientation 
stimulus distribution. We derived complete phase diagrams summarizing pattern selection in the EN model 
for fixed as well as variable retinotopy. Small interaction ranges together with low fourth moment values 
lead to either pinwheel-free orientation stripes, rhombic or hexagonal crystalline orientation map layouts as 
optimal states. Large interaction ranges together with orientation stimulus distributions with high fourth 
moment values lead to the stabilization of irregular aperiodic 0PM layouts. These solutions belong to a 
class of solutions previously called n-ECPs. This solution class encompasses a large variety of 0PM layouts 
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and has been identified as optimal solutions of abstract variational models of 0PM development [SS'. We 
showed that in the EN model due to a lack of a so-called permutation symmetry, among this family of 
solutions, states with low pinwheel densities are selected as global minima. In the extreme and previously 
unexplored parameter regime of very large effective interaction ranges and stimulus ensemble distributions 
with infinite fourth moment, permutation symmetry is restored and spatially aperiodic 0PM layouts with 
higher pinwheel density are included in the repertoire of optimal solutions. Only in this limit, the 
repertoire of optima reproduces the recently described species-insensitive 0PM design ^ and 
quantitatively matches experimentally observed orientation map layouts. None of these findings depend on 
whether the EN model is considered with variable or fixed retinotopy. 



Comparison to previous studies 

It is an important and long-standing question, whether the structure of cortical maps of variables such as 
stimulus orientation or receptive field position can be explained by a simple general principle. The concept 
of dimension reduction is a prominent candidate for such a principle (see e.g. |58(|104] for reviews) and the 
qualitative agreement between experimental data and previous numerical results from dimension reduction 



models 21,42 60,62 66 68,100,104 106 can be viewed as evidence in favor of the dimension reduction 
hypothesis. Yet comprehensive analytical investigations of dimension reduction problems and in particular 
the determination of their optimal and nearly optimal solutions have been impeded by the mathematical 
complexity of these problems. For the EN algorithm applied to the Traveling Salesman Problem, previous 
analytical results established the unselective fixed point above the first bifurcation point as well as the 



parameters at which this solution becomes unstable 107 . Subsequent work extended these results to the 
EN model for cortical map formation. The periodicity of solutions depending on the model parameters has 



been obtained by computing the eigenvalues of the Hessian matrix of the energy function 63 108 109 



Hoffsiimmer et al. confirmed these results, and computed the periodicity of the emerging patterns in the 
continuous EN model formulation by linear stability analysis of the EN gradient descent dynamics as used 
in the present study 72 . Our results extend these findings and for the first time provide analytical 
expressions for the precise layout of optimal and nearly optimal dimension-reducing maps. 
In the light of the qualitative agreement between experimental data and numerical solutions of the EN 
model previously described, it is perhaps our most surprising result that the model's optimal 
dimension-reducing maps are regular periodic crystalline structures or pinwheel-free stripe patterns in large 
regions of parameter space. In particular, the species-insensitive pinwheel statistics observed 
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experimentally 38 are not exhibited by optimal solutions of the classical EN in any of the previously 
considered parameter regimes. 

Our comparison of different numerical approaches indicates that the differences to previous studies are 
mainly attributable to differences in the sampling of the stimulus manifold in the numerical optimization 
procedures. In their seminal publication, Durbin & Mitchison used sets of 216 stimuli from the circular 
stimulus ensemble and applied a Gauss-Seidel procedure to obtain stationary configurations 21 . A similar 



procedure was used in 106 . Quite frequently, the number of stimuli used for optimization is of the same 
order of magnitude as the number of model neurons or centroids in feature space. This provides a 



relatively sparse sampling of the stimulus manifold 63 -65 . Finite stimulus sampling effects are expected 
to worsen when feature spaces of higher dimension are considered. 

The choice of small stimulus sets in previous dimension reduction studies was imposed mainly by the 
limitations of computing power. Using a parallelized implementation of the Cholesky- method for 
deterministic annealing [62]-[65l on a multicore architecture with 2TB working memory, we explored the 
dependence of the obtained near optimal solutions on the sampling of the feature space manifold over two 
orders of magnitude. We find that, the more stimuli are sampled, the closer the numerically obtained 
configurations resemble our analytical predictions. Our results on the classical EN model with 
deterministic annealing suggest that in the limit of large stimulus numbers, one would perfectly recover our 
analytical results both for periodic conditions or non-periodic boundary conditions with realistic system 
sizes. This dense stimulus sampling limit is also readily visible in our reproduction of the original Durbin & 



Mitchison sampling and the modification of the predicted map structure with stimulus number (Fig. 17). 
The finding that computational limitations prevented Durbin & Mitchison from obtaining the genuine 
predictions of their dimension reduction model should not be viewed as diminishing the importance of their 
contribution. The dimension reduction approach has played a unique and extremely productive role in 
guiding the conceptualization of cortical functional architecture. It has established an abstract view on 
cortical representations without which most of our current theoretical knowledge about candidate theories 
for cortical architectures could not have been obtained. 

Our results about optimal states of the EN for the circular and uniform stimulus ensembles however agree 
with some prior work. In 25 , the gradient descent dynamics of the EN model Eqs. ([3j|4| was used as a 
model for the emergence and refinement of cortical maps during development. Simulated visual stimulus 
features included retinotopy, orientation and eye dominance. The numerical procedures were similar to the 
one developed in the current study. Parameters were chosen such that (Is^l^) = 5.33 and a/A ^ 0.366. 



56 



This study found that an initiahy large number of pinwheels decayed via pairwise annihilation of pinwheels 
with opposite topological charge. Our analysis predicts a stripe-like OP pattern as optimal solution in this 
regime, both in the case of a fixed uniform retinotopy as well as with variable retinotopy. In our 
simulations, this state is reached after an initial phase of symmetry breaking with the generation of 
numerous pinwheels via pairwise pinwheel annihilation. Our analytical and numerical results thus confirm, 
explain, and generalize these previous findings. 

The previous results also indicated that the inclusion of eye dominance in the EN model slightly slows 
down but does not stop the pinwheel annihilation process (see 25 , Fig. 3). This raises the possibility that 
the main features of our analysis of optimal solutions for the EN model may persist when additional 
feature dimensions are taken into account. Reichl et al. in fact observed that models with interacting 0PM 
and ocular dominance maps (ODMs) exhibit a transition from pinwheel-free stripes to periodic pinwheel 
crystals similar to the transitions found in the EN 37 and demonstrated that this transition is a general 



feature of models with interacting 0PM and ODMs 110 . A rigorous characterization of map structures 



predicted by the simultaneous optimization of multiple periodic feature representations such as orientation 
preference and ocular dominance constitutes an important goal for future studies. The recent work by 
Reichl and co-workers suggests that this issue can be successfully approached using concepts from the 
nonlinear dynamics of pattern formation 37 . Finally, one recent study used the continuous formulation of 
the EN model to investigate the impact of postnatal cortical growth on the formation of ocular dominance 
columns in cat visual cortex 69 . Consistent with our results, this study also observed perfectly regular 
stripe-like patterns as stationary states in gradient descent simulations. The dynamics of the convergence 
of the ODMs towards the stripes was modified by including cortical growth into the model. However, as 
soon as growth terminated, simulated ODM layouts readily converged towards regular stripes. How 
cortical growth interacts with the formation of orientation columns is currently not understood and 
represents a further interesting topic for future studies. 

Geometric relationship between retinotopic distortions and orientation preference maps 

Experimental results on the geometric relationships between the map of visual space and the map of 
orientation preference are ambiguous. Optical imaging experiments in cat VI suggested a systematic 
covariation of inhomogeneities in the retinotopic map with singularities in the pattern of orientation 



columns in optical imaging experiments 98 . Regions of high gradient in the map of visual space 



preferentially appeared to overlap with regions of high gradient of the 0PM. In ferret, however, it has been 
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reported that high gradient regions of the map of visual space correspond to regions of low gradient in the 



0PM 67 . In tree shrew VI, no local relationships between the mapping of stimulus orientation and 
position seem to exist and the map of visual space appears to be ordered up to very fine scales 111 . In 
line with this, single unit recordings in cat area 17 revealed no correlation between receptive-field position 



scatter and orientation scatter across local cell ensembles |112 113 



Our analysis of the EN model shows that its optimal states exhibit a negative correlation between the rates 
of change of orientation preference and retinotopic position, similar to what has been observed in the 
ferret |67j. This is expected from the principle of dimension reduction and in agreement with the original 
numerical results by Durbin & Mitchison 21 . However, both in simulations of the gradient descent 
dynamics and in deterministic annealing simulations with periodic boundary conditions as well as in 
analytically obtained optimal solutions, deviations from a perfectly uniform mapping of visual space are 
surprisingly weak (see Figs. [9j [Toj [Tsj [25] in Appendix II). 

Deterministic annealing simulations with open non-periodic boundary conditions showed a substantially 
increased magnitude of retinotopic distortions. This raises the possibility that different behaviors observed 
in different experiments might be at least partially related to the infiuence of boundary effects. The 
infiuence of boundary effects is expected to decline into the interior of an area, in particular for large areas 
as VI (see [114'). In the bulk of VI, we thus expect only a weak coupling of orientation map and 
retinotopic distortions according to the EN model. In this regime, the predictions from models with 
reduced rotational symmetry (so-called Shift-Twist symmetry 115 ) about the coupling between 



retinotopic distortions and orientation preference maps 33 appear to be more promising than the weak 
effects resulting from the coverage-continuity compromise. Consistent with the measurements of Das and 
Gilbert 98 , such models predict small but significant positive correlations between the rates of change of 
orientation preference and retinotopic position 33 . Moreover, the form of the retinotopic distortions in 



such models is predicted to differ for pinwheels with positive and negative topological charge 93 . This 
interesting prediction of 0PM models with Shift-Twist symmetry deserves to be tested by measuring the 
receptive field center positions around the two types of pinwheels with single cell resolution \12\ . 

Aperiodic orientation preference maps reflect long-range intracortical suppression 

Our unbiased search through the space of stimulus ensembles with finite fourth moment revealed the 
existence of spatially aperiodic optimal solutions in the EN. It is important to realize that the selection of 
these solutions is not easily viewed as resulting from an optimal compromise between coverage and 
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Figure 20: Different patterns of evoked activity for different effective ranges of intracortical 
interaction in the EN model. The component (^^[^(x)) of the orientation map z{x.) in the direction 
of the stimulus Sz is plotted as a meshed 3D graph in a 6A x 6A patch. Color code and height of the 
projection below indicate the strength of activation. The stimulus is presented in the center of the displayed 
cortical subregion. (a) Evoked activity patterns e(x, S,z,r) for small interaction range a/A = 0.1 and 
weakly oriented stimulus with Sz = 0.01 (left) and strongly oriented stimulus with 5^ = 4 (right). rPWC 
(see upper right) are optimal in this regime, (b) Evoked activity patterns e(x, S,2:,r) for large interaction 
range a/A = 0.9 and weakly oriented stimulus with Sz = 0.01 (left) and strongly oriented stimulus with 
Sz = 4: (right). Spatially aperiodic 8-planforms (see upper right) are optimal in this regime. A uniform 
retinotopy was assumed in all cases for simplicity. 
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continuity. In fact, the continuity parameter in the respective parameter regime is so smah that solutions 



essentially maximize coverage (see Figs. 12, 15, and 16). Instead, this phenomenon reflects a different key 
factor in the stabilization of pinwheel-rich aperiodic layouts, namely the dominance of long-ranged and 
effectively suppressive interactions. This is illustrated in Fig. [20] which depicts different forms of 



stimulus-evoked activity patterns in the EN model. For a short-range interaction (Fig. 20 i), the activity 
evoked by low as well as high orientation energy stimuli is an almost Gaussian activity peak located near 
the stimulus position. The peak is shallow for low (left) and sharp for high "orientation energy" (right). In 
the corresponding parameter regime, square pinwheel crystals are the optimal solution of the EN. For a 
longer range of interaction where aperiodic 0PM layouts are the optimal states, the activity evoked by a 
single point-like stimulus is qualitatively different. Here, the activity pattern is extended and spans several 
hypercolumns (Fig. ^Ojp). It is weakly modulated for low orientation energy stimuli (left) and consists of 
several distinct peaks for high orientation energy stimuli (right). In this regime, neurons at a distance of 
several columns compete for activity through the normalization term in the EN which leads to a nonlocal 
and effectively suppressive intracortical interaction. 

It is presumably not a mere coincidence that recent studies of abstract variational models of 0PM 



development 35, 38[|95| mathematically identified this type of interaction as a key mechanism for stabilizing 
realistic 0PM layouts. It has been shown that all models for 0PM development that share the basic 
symmetries (i) translational symmetry (ii) rotational symmetry (iii) shift symmetry and (iv) permutation 
symmetry and in addition are dominated by long-range suppressive interactions, form a universality class 
that generates maps exhibiting a universal and realistic pinwheel statistics. In such models, suppressive 
long-range interactions are key to stabilizing irregular arrangement of pinwheels, which otherwise largely 
disappear or crystalized during optimization. We have stressed that the EN model as considered here obeys 
the symmetries (i)-(iii). In the limit of infinite orientation stimulus ensemble fourth moment, permutation 
symmetry (iv) is restored. The EN can thus be tuned into the above universality class by sending the 
orientation stimulus distribution fourth moment to infinity and choosing an exponentially small continuity 
parameter to realize effective long-range coupling. Indeed, the phase diagrams for abstract variational 
models of 0PM development 35 and those of the EN model found here are structurally very similar. In 
both cases, a rather large orientation stripe phase is complemented by a cascade of phase transitions 
towards more complex, aperiodic and pinwheel-rich 0PM layouts induced by long-range suppressive 
interactions. Using abstract variational models, it has been shown recently that the stabilization of regular 
crystalline pinwheel layouts can alternatively be achieved by a strong coupling between the map of 
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orientation and the map of eye dominance 37, 110 . The structure of the phase diagrams of such models 
however appears fundament ahy different from the structure of the EN phase diagrams. 
The parameter regime in which the EN model's optimal solutions exhibit the experimentally observed 
pinwheel statistics is not at all intuitive and in our opinion questions the conventional interpretation of the 
EN model to the formation of cortical feature maps. Firstly, the extremely small continuity parameter 
questions the fundamental role of a trade-off between coverage and continuity. We note that such a 
parameter regime is currently not accessible to numerical simulations. In addition, an apparently 
fundamental property for any adequate model for 0PM optimization or development, namely a 
Turing-type finite wavelength instability of the unselective state is lost in the limit r] ^ 0. At first 
sight the infinite fourth value required may appear reminiscent of the power-law distributions for 



orientation energy found in the statistics of natural images 116,117 . However, as visualized in Fig. 20 3, 
the essential property of the EN model in the infinite fourth moment regime is the occurrence of patterns 
of activity spatially extended beyond a single hypercolumn representing spatially localized point-like 
stimuli. These activity patterns mediate the long-range interactions between distant orientation columns 
which in turn cause the stability of realistic pinwheel-rich aperiodic 0PM layouts. It is obvious that 
spatially extended stimuli provide a much more plausible and realistic source of extended activity patterns 
in models for visual cortical development (for an extended discussion see [sT). Optimization models for 
cortical maps based on the representation of more complex spatially extended visual stimuli, such as 
natural scenes, rather than a model based on point-like stimuli with extreme statistics would then be a 
more appropriate basis for understanding visual cortical functional architecture. 

Comparison to the SOFM model 

Several alternatives to the EN model have been proposed as optimization approaches that can account for 
the structure of visual cortical maps. One prominent alternative dimension reduction model is the so-called 
self-organizing feature map (SOFM), originally introduced by Kohonen 59 . It is widely believed that this 
model, albeit lacking an exact energy functional ^118), implements a competition between coverage and 



continuity very similar to the EN model |56,57,66 118 . The SOFM has been reported to reproduce many 



of the experimentally observed geometric properties of visual cortical feature maps (e.g. 56,57,61 66-68]). 
The numerical procedures used in all of these studies were either the deterministic annealing procedure or 
the non-recurring application of a stimulus set without systematic assessment of pattern convergence. An 
analysis of the nontrivial stationary states of a dynamical systems formulation of the SOFM model is 
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currently lacking. The main difference between the SOFM model and the EN model is that the activation 
function by definition has the form of a stereotypical Gaussian and competition is incorporated by a hard 
winner-takes-all mechanism. As a consequence, it is not obvious that a long-range suppressive interaction 
regime can be realized in this model. According to our analysis, one would thus expect orientation stripes 
and rhombic pinwheel crystals as nontrivial stationary states of the SOFM model. In a very recent study 
of the SOFM algorithm that used a numerical procedure similar to the gradient descent simulations 
developed in this paper, both pinwheel annihilation and rhombic pinwheel crystallization have been 
observed [119]. In addition, one study that examined the SOFM model for orientation and retinotopy 
found a fast convergence to pinwheel-free stripe-like solutions for a wide parameter range [25] . In view of 
these results, it seems worthwhile to also reexamine the SOFM model with respect to its stationary states. 

Rugged or Smooth Energy landscape 

As for many optimization problems in biology, the optimization of visual cortical functional architecture 
has been considered a problem characterized by a rugged energy landscape 77 . In case of the EN model 
the expectation of a rugged energy landscape at first sight seems quite plausible. Originally, the Elastic 
Network algorithm was invented as a fast analogue method to approximately solve NP-hard problems in 
combinatorial optimization such as the traveling salesman problem (TSP) 78j[79j. In the TSP, the stimulus 
positions correspond to the locations of cities a salesman has to visit on the shortest possible tour. In 
problems such as the TSP, the energy functionals to be minimized are known to possess many local minima 
and the global minimization of these functionals generally represents an extremely difficult problem \79\ . 
Our analysis reveals that the trade-off between coverage and continuity for the mapping of a continuous 
feature space manifold leads to a much simpler structure of the energy landscape. This is also indicated by 
the fact that almost all of our gradient descent dynamics simulations readily converged to the predicted 
global minimum of the energy functional. Fig. [2l] illustrates the smooth structure of the EN energy 
landscape close to pattern formation threshold for different model parameters for a one dimensional path 
through the state space. In this landscape, the small set of stable planforms correspond to local minima of 
the EN energy functional, and unstable planforms to saddle points in the energy landscape. The optimal 
states correspond to global minima. Note that along the depicted state space path, unstable stationary 
solutions may appear as local minima if the unstable directions along which the energy decreases are 
orthogonal to the path. 

What is the origin of this qualitative difference in the shape of the energy landscapes? In the traveling 
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Figure 21: Illustration of the EN energy landscape close to pattern formation threshold. The 

variation of the energy between states of ideal orientation stripes (OS), the 2-ECP state, square pinwheel 
crystals (sPWC) and possible mode configurations for 8-ECPs is shown for the case that (a) the OS state 
has the lowest energy (54 = 0, a/A = 0.3), (b) the sPWC state has the lowest energy (54 = 0, a/ A = 0.1), 
and (c) the most anisotropic 8-ECP has lowest energy (54 = 100, <j/A = 0.8). The energy values between 
the state are computed from a state obtained by linear interpolation between two neighboring states on the 
X-axis. Note that not all local minima in a-c correspond to a stable fixed point of the amplitude dynamics 
(see text). 
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salesman problem, the finite repertoire of possible tours consists of all permutations of the N cities that the 
salesman has to visit. Via self-organized competition between the aim to visit all cities and the aim to 
minimize the path length, the elastic network algorithm converges to a specific ordering of the cities that 
eventually yields a very short tour. Most likely, the qualitative difference to the EN model for visual 
cortical map architecture originates from the transition from a finite number of cities to a continuum. 
When the elastic network algorithm is considered with an ensemble of cities (or stimuli) distributed 
according to a continuous probability density function, there is no discrete repertoire of tours. Both, the 
repertoire of tours as well as the path through the landscape of cities or equivalently the space of visual 
stimulus features are determined by self-organization. The first is generated by the symmetry breaking 
mechanism that leads to the instability of the homogeneous state. The second corresponds to the selection 
of one of the many nontrivial stable steady states. 

An interesting property of the EN model dynamics that can be inferred from the energy landscape 



depicted in Fig. [21] is the type of competition between two stable stationary states, where both are present 
in the system with a wall or a domain boundary between them. The motion of the wall or domain 
boundary is predicted to proceed in the direction that increases the fraction of the pattern with lower 
energy. An example of such competition can be seen in Fig. [23^ . At t = lOOr, a small domain with an 
sPWC state is present. The area of this region is gradually reduced over the time course of the simulation 
until the pinwheel-free optimal state is reached. 

Are simple 0PM layouts an artifact of model simplicity? 

The perfectly periodic types of stationary solutions (stripes, crystals) that appear to dominate the classical 
EN model for retinotopy and orientation have been found in other models of visual cortical layouts that are 
relatively abstract. One might therefore suspect that they represent a mere artifact of model simplicity. 
One conceptually appealing approach where perfectly periodic layouts have been found is wiring-length 
minimization 27 . According to this hypothesis, the structure of an 0PM can be understood by 
minimizing the total length of dendritic and axonal processes. Maps obtained by stopping minimization of 
wire length exhibited qualitatively realistic layouts (see Fig. 6 in 27 ). Complete optimization, however, 
leads to either stripe-like pinwheel-free patterns or rhombic pinwheel crystals, identical to the ones 
obtained in our investigation of optimal solutions of the EN model 27 . Similarly, stripe-like and rhombic 



optima have been found in several abstract vector-field approaches for 0PM development 31 , 33 120 



It is ruled out by two observations, that the crystalline and perfectly periodic optima observed in all four 



64 



optimization models, the EN model, the SOFM, the wiring-length minimization model, and the vector-field 
models are a mere artifact of the abstract order parameter field description of cortical selectivity patterns 
that is common to these approaches. Firstly, equally simplistic order parameter models for 0PM 
development with long-range interactions have been shown to reproduce spatially irregular map 



layouts [35 38,95 . The occurrence of periodic optimal solutions is thus not a necessity in this model class. 
Secondly, pinwheel crystallization has also been observed in detailed network models for the development 
of orientation preference maps, notably in the first ever model for the self-organization of orientation 



selectivity by von der Malsburg in 1973 14, 121 . Thus, on the one hand the phenomenon of pinwheel 
crystallization is thus not restricted to simple order parameter models and on the other hand abstract and 
mathematically relatively simple models can exhibit complex and biologically realistic optimal solutions. 

Map rearrangement and layout optimization 

Irrespective of the optimization principle invoked to describe the structure of visual cortical maps, several 
common features of the resulting dynamics have been observed. The dynamics of optimization models 
usually starts with a phase of pattern emergence, where selectivity to visual features arises from an initially 
homogeneous unselective or weakly selective state. As we and others have shown, feature maps in these 
models continue to evolve after single cell selectivities reach mature levels. In fact, the phase of initial 
pattern emergence is typically followed by a prolonged phase of rearrangement of selectivities and 
preferences until a stable configuration is reached that represents a genuine optimum. This is not an 
exceptional type of dynamics but rather constitutes the generic expectation for a spatially extended 
system 



84,85 



What drives the second phase of map rearrangement? The initial emergence of feature selectivity is 
predominantly a local process in which merely neighboring units interact with each other to roughly match 
their selectivities. In the resulting spatial layout, selectivities are therefore far from being optimally 
arranged in space with respect to the global organization of selectivities on larger scales. Depending on the 
interactions incorporated in the model, local matching processes may (i) effectively propagate through 
space optimizing the pattern over gradually increasing spatial scales or (ii) distant sites may start to 
directly interact with each other to guide a rearrangement towards a globally optimized pattern after their 
initial emergence. 

An illustrative example is provided by the emergence of pinwheel- sparse orientation stripes. Qualitatively, 
it is easy to see that a pattern of orientation stripes satisfies the continuity constraint very well. In a stripe 
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pattern, preferred orientations are constant along one direction in space, realizing the absolute minimum of 
the orientation gradient in this direction. Reaching such a configuration obviously requires to select the 
preferred orientation at widely separated sites (along the stripe axis) to be identical. Because initially such 
sites develop independent preferred orientations, the optimized column layout can only emerge through a 
secondary rearrangement process. If the dominant low energy state has low pinwheel density, the later 
phase is governed by pinwheel motion and pairwise pinwheel annihilation. If this state is pinwheel-rich, e.g. 
a pinwheel crystal or an aperiodic pinwheel-rich state, both pinwheel annihilation and pinwheel creation 
together with a coordinated rearrangement of pinwheels are expected to occur. 

The local, essentially random processes during the initial emergence of a first pattern are in principle 
incapable of directly generating an optimized layout. In fact, it has been established that this initial 
so-called symmetry breaking phase will in general produce a random arrangement of selectivities of 
model-insensitive statistics ^ 32 , 36| . The occurrence of some form of secondary reorganization is thus a 
qualitative prediction of any optimization model, provided that the optimal map is not seeded by an innate 
mechanism. The results presented in this study and many reports demonstrate that Hebbian plasticity is 
capable and often expected to achieve such rearrangements. 

In gradient descent dynamics simulations of the EN model for retinotopy and stimulus orientation with 
conventional stimulus ensembles, pinwheel densities were found to be strongly time-dependent after the 
initial column formation (see e.g. Figs. [6j[7|). In particular, the timescale for the establishment of full 
orientation selectivity and the time needed for either annihilation of a substantial fraction of pinwheels or 
their crystallization into periodic pinwheel crystals are in the same range of tens of tau. A similar time 
dependence of pinwheel density has also been observed in other models for 0PM development with 



periodic optima |i35,38 . Pinwheel annihilation in the EN can be slightly slowed down by additional 
features such as retinotopy Figs. 10, 17 or ocular dominance 25l[37l but not by orders of magnitude. For 
this reason, signatures of the periodic optima of a developmental dynamics become visible at rather early 
simulation stages. Long-term minimization is apparently not essential to express the main layout features 
of the global minimum. 

Because the main features of the dominant optimal solutions become apparent immediately after 
orientation selectivity saturates it appears not easy to reproduce the species-independent map layout in 
models with periodic crystalline optima by pattern freezing. In our simulations to match even only the 
pinwheel density, a very precise timing of the freezing point would be required. There is currently no 
evidence for such a freezing mechanism in early development. In cats and ferrets, cortical maps for ocular 
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dominance, orientation or direction arise on a timescale between hours and a few days (e.g. 13, 122, 123 ). 
The underlying circuits can be rapidly modified, e.g. by deprivation experiments, even on the timescale of 
hours [124[[125| weeks after full selectivity has been established. Recently, evidence for long-term visual 
cortical circuit reorganization after the emergence of feature selectivity during normal development has 
emerged in diverse systems. In mouse, for example, activity-dependent changes induced by normal visual 
experience during the critical period, i.e. long after the primary emergence of orientation selectivity, have 



been shown to gradually match eye-specific inputs in the cortex 126 . Specifically, the data from mouse 
indicates that preferred orientations in the two eyes initially often emerge unmatched and subsequently 
change towards one binocularly matched orientation preference. Because preferred orientations in the two 
eyes initially are statistically independent, this suggests that neurons can rotate their orientation 
preferences up to at least 45° during postnatal development. This is reminiscent of pairing experiments in 
kitten visual cortex in which Fregnac and coworkers induced neurons to changed their preferred orientation 
by up to 90° after pairing of a visual stimulus with intracortical stimulation 127[|128| (see also [129] ). Also 
in the cat, visual cortical orientation columns in visual areas VI and V2 have been found to undergo 
rearrangement during the late phase of the critical period 41 . In this process, columns in mutually 
connected regions of areas VI and V2 or in retinotopically matched regions in left and right hemisphere 
areas become progressively better matched in size. In the same species, a systematic reorganization of 
ocular dominance columns during postnatal development has been observed [69] . Essential features of this 
columnar rearrangement are reproduced by the EN model for ocular dominance patterns simulated in a 
growing domain. 

In view of these observations, it seems unlikely that aperiodic orientation maps in the visual cortex 
represent frozen transient states of a developmental dynamics whose attracting layouts are pinwheel 
crystals or pinwheel free states. In fact, models for the activity-dependent development of OPMs with 
aperiodic optima predict only subtle changes of the 0PM layout during the convergence after the 



establishment of selectivity 35,38 . This might also explain the apparent stability of cortical maps during 



normal development over short periods 122 . Further studies of the long-term rearrangement and 
stabilization of cortical functional architecture are needed to exhaustively characterize such processes. 
Given the fundamental role of map reorganization for any optimization theory of visual cortical 
development, chronic imaging experiments tracking the spatial arrangement of feature selectivities in 
individual animals beyond the emergence of selectivity and through later developmental stages are 
expected to be highly informative about fundamental principles of visual cortical optimization. 
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Conclusions 

Together with recent progress on the quantitative characterization of cortical functional 
architecture [38l[69{[95] , the current study lays the foundation for a mathematically rigorous and biologically 
informative search for optimization principles that successfully explain the architecture of columnar 
contour representations in the primary visual cortex. A mathematically controlled and quantitatively 
precise determination of the predictions of candidate optimization principles is demanded by accumulating 
evidence indicating that geometrical features of visual cortical representations are biologically laid down 



with a precision in the range of a few percent 38 130 131 . Such data is expected to substantially reduce 
the range of candidate optimization principles that are consistent with biological observations. In 
particular, for the principle that cortical orientation maps are designed to optimally compromise stimulus 
coverage and feature continuity, our analysis demonstrates that the classical EN model for orientation 
preference and retinotopy essentially fails at explaining the biologically observed architecture. Our finding 
that the EN model exhibits biologically realistic optima only in a limit in which point-like stimuli are 
represented by complex spatially extended activity patterns corroborates that large-scale interactions are 
essential for the stabilization of 0PM layouts with realistic geometry ^35j[39j[88j[95] . In the light of these 
results, principles for the optimal representation of entire visual scenes by extended cortical activity 
patterns appear as promising candidates for future studies (see also [54j). In fact, there is recent evidence 
that visual cortical activity becomes progressively better matched to the statistics of natural stimuli but 
not to simplistic artificial stimulus ensembles 132 . We expect the methods developed here to facilitate a 
comprehensive characterization of such candidate principles. 
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Methods 

Expansion of Elastic Net equation 

In order to analytically calculate the approximate optimal dimension-reducing mappings in the EN model 
with fixed retinotopy, an expansion of the nonlinear EN 0PM dynamics Eq. (|3| up to third order around 
the unselective fixed point has to be derived. This expansion is briefly sketched in the following. 
Eq. (|3| with r(x) = is of the form 

dtz{yi, t) = M^[z] + r]Az{yi, t) , 

where My^z] is a nonlinear functional of z{-^t)^ parametrized by the position x. Clearly, the diffusion term 
contains no nonlinear terms in z{-^t) and therefore third order terms of the dynamics 9tz(x, t) exclusively 
stem from third order terms of the Volterra series expansion of the functional A^x [z] around the fixed point 
z(x, = 0. By the symmetry Eq. ([s]), only third order contributions of the form N^[z^ z] are allowed, i.e. 
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Collecting all the terms yields 
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and 

i^4(xi,X2,X3,X4,X5,X6) = (^'+^2+x^+x^+x^+x^)/(8^') . 

The coefficients aj for various orientation stimulus ensembles are given in the Results section. 
Adiabatic Elimination of r(x, t) 

In order to analytically calculate the approximate optimal dimension-reducing mappings in the EN model 
with variable retinotopy, an expansion of the nonlinear EN retinotopy and orientation map dynamics Eqs. 
(|3j [4| up to third order around the nonselective fixed point has to be derived and retinotopic distortions 
have to be adiabatically eliminated. Both of these calculations are briefly sketched in the following. Eq. 
(|3| is of the form 

dtz{^, t) = Af^ [z, r] + 7?Az(x, t) , 

where JVx[z^r] is a nonlinear functional of z{-^t) and r(-,t), parametrized by the position x. The diffusion 
term contains no nonlinear terms in z{-^t) and therefore third order terms of the dynamics of ^(x, t) 
exclusively stem from third order terms of the Volterra series expansion of the functional A/'x[^,r] around 
the fixed point {z{'x,t) = 0, r(x, t) = 0}. By the symmetry Eq. ([8|, only terms in form of a cubic operator 
Ns[z^ z] and a quadratic operator Q^[r^ z] are allowed when expanding up to third order. Ns[z^ z] is 
given in Eq. ( [37[ ). (5[^,r] can be calculated via 

and this yields 

/ '^''y ^(y)K2(y - x)) + j d'y (r(y),^(y)K^(y - x)) 

^ d^yd^w;2;(y)(r(w),K3(y-x,w-x,y-w)) , 
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where (•, •) denotes the scalar product between two vectors and 

K^(x) = e-"'/4^'x (39) 
K5(xi,X2,X3) = e- ' e/- ' [xi + X3] . (40) 
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In complete analogy, by expanding the right hand side of the dynamical equation for the retinotopic 
distortions (Eq. Q) up to second order, the vector- valued quadratic operator Q^[z,z] can be obtained as 
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Inserting r(x) = — -^[Q^[z, z]] into (5^[r, z] and using the linearity of ^ as well as the bilinearity of 
both, Q^[z, z] and Q^[r, z], yields a sum N^[z^ = Xljii a^N^^ with 
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The coefficients are given by 
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where the second equal sign is vahd for (l^^p) = 2. 
Amplitude Equations from N§[z,z,z] 



We catalog the numerous stationary solutions of Eq. (13) following 35 , by considering planforms 

N 



with an even number N of modes with amplitudes Aj and kj = A:c(cos(27rj/A^), sin(27rj7A^)). In the 
vicinity of a finite wavelength instability - where the nonselective state z(x) = becomes unstable with 
respect to a band of Fourier modes around a finite wave number kc - by symmetry, the dynamics of the 
amplitudes Aj at threshold has the form 

AT AT 
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where j denotes the index of the mode antiparahel to the mode j, kj = — k^-, and the coefficients 
gij = (1 - ^Sij)g{\ai - aj\) and fij = (1 - % - Si-j)f{\ai - aj\) only depend on the angle \ai - aj\ 
between mode i and j. The angle-dependent interaction functions g{a) and /(a) are obtained from Eq. 

as 



( 13 ) by a multi scale expansion 35 84 85 
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(43) 
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where ko = /cc(l,0) and h(a) = /Cc(cosa, sina). /(a) is 7r-periodic, since the right hand side of Eq. (44) is 
invariant with respect to the transformation h(a) h(a + tt) = — h(a). g{a) is 27r-periodic in general. If, 



however, the nonlinearity is permutation symmetric (Eq. ([34])) it can be seen from Eq. (43) that g{a) is 
TT-periodic as well. 



Stability of stationary planform solutions 



To determine the stability of fixed points of the amplitude equations Eq. (42), the eigenvalues of their 



stability matrices have to be determined. In general, for any fixed point A = of the dynamical system 
A = F(A) with complex- valued A and F, we have to compute the eigenvalues of the Hermitian 2N x 2N 
matrix 

/ dF dF 
V dA dA 

For the system of amplitude equations, we obtain 



A=Ao 



OF, 



N 



rSik-Sik \ ^gij\Aj\'^ - AiQikAk - Ai- fik{Af,- ^ Ak) 



N 



AigikA}^ ^i-fc I ^ ^ fij^j^j 



Stability of a solution, or more precise intrinsic stability is given, if all eigenvalues of M are negative 
definite. Extrinsic stability is given, if the growth of additional Fourier modes is suppressed. To test 
whether a planform solution is extrinsically stable, we introduce a test mode T such that 

N 
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with k/3 = (cos/3, sin/3)/Cc. We insert this ansatz into Eq. (15) and obtain 



N 



3 

as the dynamics of the test mode T, where g{l3) is the angle-dependent interaction function corresponding 
to N^[z^ z]. For the solution T = to be stable, we therefore obtain the condition 

N 

r-Y.9{p- Pj)\Aj\^ < 0, Va G [0,2^] , 

3 

where we assumed ^ kj,k~. These conditions for intrinsic and extrinsic stability were numerically 
evaluated to study the stability of n-ECPs and rPWCs. 

Coupled Essentially Complex Planforms 

In the Results section, we presented a closed form expression for the retinotopic distortions associated via 



Eq. (28) with stationary planform solutions of Eq. (29). Here, we sketch how to explicitly calculate this 



representation. We start with the ansatz 

n 



(45) 



for the orientation preference map z(^). Note that this general ansatz includes essentially complex 
planforms as well as rhombic pinwheel crystals. To simplify notation, we denote the individual terms in 



Eq. (41) 



Wna^ ^^^^ 

d2yKS(y-x)|z(y)|2 



327rc76 



// 



d yd i(;K3(y - x,y - w, w - x)z(w)z(y) . 



Each of the Q^^, ^ = 1, 2, 3 can be evaluated for the ansatz Eq. (45) and we obtain 

{\s^?)e-K ~ 



Qi 



-ikfcx 



2^2 E (^i^O (ki - kfe) sin ((k, - kfc) x) 



{AjAk) (k^- - k/e)cos((k^- - k/e)x)} 



Q2 



Zee 



{2a^-{\s.?)) " 

j<k 



2a2 



^ g-(k,-k,) (k, -kfc) {A,Au) sin ((k, - kfc) x) 



+^ {AjAk) cos ((k^- - kfc) x)} 
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Q3 



^^^^ E(^^- - ^k)e-^^ {J? (A, A,) sin((k, - k,)x) 



j<k 



(A^ Afc) COS ((k^- - k/e) x)} . 



All resulting terms are proportional to either (k^ — kj) sin((k^ — kj)x) or (k^ — kj) cos((ki — kj)x) , i ^ j. 
These functions are longitudinal modes (see Fig. ^) which have been identified as eigenfunctions of the 
linearized dynamics of retinotopic distortions L^[r] with eigenvalue 

A2(|k, - k,|) = -|ki - k^firjr + e-'^'l'^'-'^^lV^) . 
Hence, they are eigenfunctions of L~^[r] with eigenvalue 1/A£(|ki — kj|). Using this when inserting in Eq. 



(28) and setting (Is^P) = 2, we obtain expression (31) for the retinotopic distortions belonging to an 
arbitrary planform. 



Phase diagrams 

To compute the regions of minimal energy shown in Figs. [6| p^Oj [23| [25| (Appendix) , 12, 15, and 16, we first 



computed the fixed points of Eq. (42) at each point in parameter space. For n-ECPs, we constructed the 



coupling matrix g in Eq. ( 22 ) for all mode configurations not related by any combination of the symmetry 



operations: (i) Translation: Aj A^e^^^^, (ii) Rotation: Aj A^+aj, (hi) Parity: Aj A(7v-j)-- Via 



Eq. (22), we then computed the absolute values of the corresponding amplitudes. If Yl^^i (s ^) ij ^ for 



all z, a valid n-ECP fixed point of Eq. (|42| was identified. Its energy was then computed via Eq. ( |23[ ). For 
orientation stripes and rhombic pinwheel crystals, the derived analytical expressions for their energy Eqs. 



(18, 20) were evaluated. To analyze the stability of the fixed points, the conditions for intrinsic and 
extrinsic stability (see above) were numerically evaluated. 



Numerical Procedures - Gradient descent simulations 

To test our analytical calculations and explore their range of validity, we simulated Eqs. ([3| and (|4| on a 
64x64 grid with periodic boundary conditions. Simulated systems were spatially discretized with typically 
8 grid points per expected column spacing Amax of the orientation preference pattern (see Results section) 
to achieve sufficient resolution in space. Test simulations with finer discretization (16 and 32 grid points 
per Amax) did not lead to different results. Progression of time was measured in units of the intrinsic 
timescale r (see Results section) of the pattern formation process. The integration time step 5t was 
bounded by the relevant decay time constant of the Laplacian in Eq. (|3| around kc and by the intrinsic 
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timescale r of the system, using 5t = min {l/(207^A:^), r/lO}. This ensured good approximation to the 
temporahy continuous changes of the patterns. We used an Adams-Bashforth scheme for the first terms on 
the respective r.h.s. of Eqs. (3j4). The second terms (diffusion) were treated by spectral integration, 
exhibiting unconditional numerical stability. The stimulus positions were chosen to be uniformly 
distributed in retinal coordinates. The stimulus averages in Eqs. (|3j |4| were approximated by choosing a 
random representative sample of Ns stimuli at each integration time step, with 

where n corresponds to the dimensions of the feature space in addition to the two retinal positions (in our 
case, n = 2), = (L/Amax)^ the squared aspect ratio of the simulated system in units of A^, Eg the 
resolution in feature space, TVq the number of stimuli we required to approximate the cumulative effect of 
the ensemble of stimuli within each feature space voxel With A^o = 100 and £s = 0.05, we ensured a 

high signal to noise ratio for all the simulations. Typical values for Ns were between 2.5x10^ and 4x10^. 
All simulations were initialized with z(x, t = 0) = 10~^e*^^^^^^ and r(x, t = 0) = 0, where the ^(x) are 
independent identically distributed random numbers uniformly in [0, 1]. Different realizations were 
obtained by using different stimulus samples. 

Stimuli were drawn from different distributions, each with (I^^P) = 2. We considered (i) stimuli uniformly 
distributed on a ring with \sz\'^ = a/2 (circular stimulus ensemble) (ii) stimuli uniformly distributed within 
a circle {s^, l^^l < 2} (uniform stimulus ensemble) and (iii) a Gaussian stimuli ensemble with 
Ps^ = l/(27r) exp(— |52p/2). In addition, we considered mixtures of a high-fourth moment Pearson type VII 
distribution and the circular stimulus ensemble. The Pearson distribution is given by 

1 



Psz — _ r^/.__ 1 In 
2 ' 2^ 



aB{m 



where 5(-, •) is the Beta function 133 and a = \/2m — 3, and m = | + such that (l^^p) = 2, 
(Is;^!'^) = 7 or equivalently 54 = 7 — 4. 

In addition to simulations in which independent sets of stimuli were used for evaluating the stimulus 
average in Eqs. ([sj [4]) for every time step, we also performed simulations in which the same fixed set of N 
stimuli was used (see Results). To determine the time step 5t in these simulations, we first calculated 



(parameters as in regular simulations) which yields the number of stimuli presented to the model in one 
intrinsic time unit r in regular simulations. To subject the network to the same number of stimuli per 
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intrinsic time scale r in fixed stimulus set simulations, the integration time step 6t was in this case chosen 
as 

For small this resulted in very small integration steps. For very large time steps were identical to 
the regular simulations. Different realizations were obtained by different but fixed stimulus sets. In all 
simulations with fixed stimulus sets, stimuli were drawn from the circular stimulus ensemble. All other 
numerical methods were chosen as in regular simulations. 



Numerical Procedures - Solving the EN model with deterministic annealing 

A large body of previous work has solved the EN models for various aspects of visual cortical architecture 
for discrete fixed sets of stimuli and using deterministic annealing. We therefore also used deterministic 
annealing with fixed discrete sets of stimuli to solve the EN model for the most frequently used stimulus 
distribution. This allowed us to better compare our analytical and numerical results based on the gradient 
descent dynamics for a continuous stimulus with prior results. For the discrete deterministic annealing 
approach, cortical maps are described by a collection of M centroids {ym}m=i ^^^^ 
represented as a I) x M matrix Y = (yi, . . . , Ym)- Maps forming a compromise of coverage and continuity 
are obtained for {x^}^^;^ represented diS a D x N matrix X = (xi, . . . , xm). In our case, d = 4. The 

trade-off between coverage and continuity is formulated by the energy function 

E{Y,a) = -aa^log ^ g-^ll'^'^ir + ^tr (Y^YS) . (46) 



The matrix S determines the topology of the network as well as the boundary conditions and is typically 
derived from a discretized derivative based on a finite difference scheme or stencil (see below). For large N 



and M, the energy function Eq. (46) is equivalent to the energy functional of the continuum formulation 
Eq. ([T]) for /3 = 77 and S implementing the discretized Laplacian operator in two dimensions. 
Following 62-65,71 we minimized the EN energy function Eq. ([T]) by an iterative deterministic annealing 
algorithm, starting with a minimization for large a and tracking this minimum to a small value of a. As in 
[4], we reduced a from ainit = 0.2 to the point at which the amplitude of the orientation maps saturate 
(cr ^ 0.03), following a = ainit x where j counts the annealing step. This choice tracks stationary 
solutions of the EN to parameters that are very far from threshold. For high precision tracking of 
solutions, we used an annealing rate of x = 0.999. 
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At each value of a, setting the gradient of Eq. (46) to zero yields a nonlinear system of equations 



YA = XW with A = G + a/3 ( ) , (47) 



2 



Here, the x M-matrix W is given by 



_ 1 II -y 

e 2 II 



and the M x M-matrix G is 



= 



N 

n=l 

A is a symmetric positive-definite M x M matrix. The M x M matrix A is symmetric and 



positive-definite. Since both, G and W depend on Y, this equation is nonlinear in Y and has to be solved 



iteratively. Following 62 - 65 71 , we solved eq. (47) at each value of a and for each iteration directly via 
Cholesky-factorization. 

We implemented periodic and non-periodic boundary conditions by appropriate choice of the matrix S. S 
must be positive (semi) definite for the energy to be bounded from below. We used the canonical 2D 
Laplacian stencil of order 2, to construct the M x M matrix 

/-4 + 2a 1 ••• 1 OO--- 1 \ 

1 -4 + a 1 ••• 100 ••• 

1 -4 + al ... 1... 

1 -4 2a 1 

1 -4 + a 1 

••• 1 -4 1 ••• 1 

1 -4 1 ••• 1 

1 -4 + a 1 

1 -4 + a 1 
V 1 1 -4 + 2a / 

Here, a = for periodic boundary conditions and a = 1 for non-periodic boundary conditions. In Appendix 
II, we also present simulation results for a 4th derivative stencil, in which S'^ was used for the continuity 
term. We used random stimulus positions and orientations as well as stimuli arranged on a grid in feature 
space. For random stimuli, positions were drawn from a uniform distribution in [0, 1] x [0, 1]. Orientations 



Sz were drawn from the circular stimulus ensemble, with l^^l = 0.08 as in 65 . Stimuli from grid-like 
ensembles were distributed evenly-spaced in [0, 1] x [0, 1] and contained 2^ evenly space orientations with 
IsJ =0.08. 



78 



To compute the energy of pinwheel-free configuration, we initiated simulations with a stripe-like 
orientation preference pattern with the same typical spacing as the observed orientation maps and 
annealed from a = 0.035 to a = 0.03. 

To enable comparison between simulations with different numbers of stimuli, we scaled the continuity 
parameter proportionally to N such that the equivalent r] was approximately constant. The simulated 
domain then contained a comparable number of hyper columns for all stimulus numbers. 



Pinwheel density from simulations 

Pinwheels locations in models were identified by the crossings of the zero contour lines of real and 
imaginary parts of the orientation map. Estimation of local column spacing A(x) was done using the 



wavelet analysis introduced in 130,131 . In short, an overcomplete basis of complex Morlet wavelets at 
various scales and orientations was compared to the 0PM pattern at each location. A(x) was estimated by 
the scale of the best matching wavelet. The mean column spacing (A(x))^ of a given map was then 
calculated from the local column spacing by spatial averaging. For details we refer to [38lpQl[T3l] . Given 



Np^ pinwheels in a simulated cortical area of size L , we defined the pinwheel density ^ 38,104 



p = Np. 



The pinwheel density p is a dimensionless quantity and depends only on the layout of orientation columns. 
The pinwheel density defined in this way is large for patchy and small for more band-like columnar layouts. 
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Appendix I 

The impact of non-oriented stimuli 

The main text of the present paper contains a complete analysis of optimal dimension-reducing mappings 
of the EN model with a circular ensemble of orientation stimuli. These optima are simple regular 
orientation stripes or square pinwheel crystals. The circular orientation stimulus ensemble, however, 
contains only stimuli with a fixed and finite "orientation energy" or elongation \sz\' This raises the 
question of whether the simple nature of the circular stimulus ensemble might restrain the realm of 
complex dynamics in the EN model. The EN dynamics are expected to depend on the characteristics of 
the activity patterns evoked by the stimuli and these will be more diverse and complex with ensembles 
containing a greater diversity of stimuli. Therefore, we also examined the EN model in detail for a richer 
ensemble of stimuli. In this ensemble, called a uniform stimulus ensemble in the following, orientation 
stimuli are uniformly distributed on the disk {s^^ \sz\ < 2}, a choice common to many previous studies. 



e.g. |19 25 , 82||94| . The uniform ensemble in particular contains unoriented stimuli with \sz\ =0. 
Intuitively, the presence of these unoriented stimuli might be expected to fundamentally change the 
importance of pinwheels in the optimal 0PM layout. Pinwheels' population activity is untuned for 
orientation. Pinwheel centers may therefore acquire a key role for the representation of unoriented stimuli. 
As such an effect should be independent of retinotopic distortions and to aid comparison with our previous 
results, we will again start with a fixed uniform retinotopy r(x) = 0. 

The linear stability properties of the unselective fixed point are independent of the ensemble of orientation 
stimuli ((l^^p) = 2 throughout this paper). The coefficients in Eq. ([l4|, however, of course depend on the 
fourth moment of the stimulus distribution. Inserting (l^^l^) = 16/3 into Eq. 



32 



we obtain 
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Figure 22: Angle-dependent interaction functions for the EN model with fixed retinotopy and 
uniform orientation stimulus ensemble, (a, b) g{a) and f{a) for a/ A — 0.1 (a) and <j/A = 0.2 (b). 
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Figure 23: Optimal solutions of the EN model for uniform stimulus ensemble and fixed repre- 
sentation of visual space, (a) At criticality, the phase space of this model is parameterized by either the 
continuity parameter r] (blue labels) or the effective interaction range a/ A (red labels, see text), (b-c) OPMs 
(b) and their power spectra (c) in a simulation of Eq. (|3| with r(x) =0, r = 0.1, a/A = 0.12 (r^ = 0.57) and 
uniform stimulus ensemble, (d) Analytically predicted optimum for a/A < 0.15 (rhombic pinwheel crystal), 
(e) Pinwheel density time courses for four different simulations (parameters as in b; gray traces, individual 
realizations; black trace, simulation in b; red trace, mean value), (f) Mean squared amplitude of the station- 
ary pattern in simulations (parameters as in b) for different values of the control parameter r (black circles) 
and analytically predicted value (solid green line), (g-h) OPMs (g) and their power spectra (h) obtained in 
a simulation of Eq. ([3| with r(x) = 0, r = 0.1, a/A = 0.15 (77 = 0.41) and uniform stimulus ensemble, (i) 
Analytically predicted optimum for a/ A > 0.15 (orientation stripes), (j) Pinwheel density time courses for 
four different simulations (parameters as in g; gray traces, individual realizations; black trace, simulation in 
g; red trace, mean value), (k) Mean squared amplitude of the stationary pattern in simulations (parameters 
as in g) for different values of the control parameter r (black circles) and analytically predicted value (solid 
green line). 
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The angle-dependent interaction functions are then given by 

+ _!_ (^^2klaHoosa-l) _ ^ ^e"'"'"' sinh^ {l / 2klcj^ COS a) 

f{a) = ^ (l - {cosh{2kla\osa) + 2cosh{k^^a\osa)) + 2e-''>^^ 

+ ^ (e-2'=>' cosh(2fc2a2 cos a) - l) + ^e-^>'>" sinh^ (l/2fce'^2 cos a) . 

Both functions are depicted in Fig. [22] for two different values of the effective intracortical interaction 
range a/ A. They qualitatively resemble the functions depicted in Fig. [s] Fig. 23 displays the phase 



diagram of the EN model with uniform stimulus ensemble. As summarized in the main part of the paper, 
it is almost identical to that obtained for the circular stimulus ensemble (Fig. |6|. Two different optimal 
states are found, square pinwheel crystals (sPWC) and orientation stripes (OS) separated by a phase 
border at a/A ~ 0.15. Both fixed points are stable for all a/A. Fig. [23|3-k demonstrates, that these 
analytical results are confirmed by direct numerical simulations of Eq. ([3| with r(x) = 0. As for the 
circular stimulus ensemble, we also tested the stability of stationary n-ECP solutions with 2 < n < 20 by 
numerical evaluation of the criteria for intrinsic and extrinsic stability (Methods). We find all n-ECPs with 
2 < n < 20 intrinsically unstable for all interaction ranges a/A. The simple phase space structure 
furthermore apparently remains unchanged if we consider the model far from pattern formation threshold 
as shown in Fig. [24j Simulations bear a close resemblance to the simulations with circular orientation 
stimulus ensemble (Fig. [t]). Either convergence to sPWC-like patterns or patterns with large orientation 
stripe domains is observed. Again, pinwheel annihilation in the case of large a/ A is less rapid than close to 
threshold (Fig. [24^,b). The linear pinwheel- free zones increase their size over the time course of the 
simulations, eventually leading to a stripe pattern. For smaller interaction ranges a/A, the 0PM layout 
rapidly converges towards a crystal-like rhombic arrangement of pinwheels with dislocations and pinwheel 
density close to 4. 

Fig. [25] shows that taking retinotopic distortions into account yields an almost identical picture compared 
to the circular stimulus ensemble. For small interaction range a/ A, the analytically predicted optimum is a 
quadratic pinwheel crystal with pinwheel density p = 4. For larger a/ A, the analytically predicted 
optimum is an orientation stripe pattern with pinwheel density p = 0. Our results are confirmed by direct 
simulations of Eq. ([3j[4| (Fig. [25|3-e). The simulation results are virtually indistinguishable from the 
circular stimulus ensemble. 
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Figure 24: Numerical analysis of the EN dynamics with uniform orientation stimulus ensemble 
and fixed representation of visual space far from pattern formation threshold, (a) OPMs and 
their power spectra in a representative simulation of Eq. (|3| with r(x) = 0, r = 0.8, a/A = 0.3 (r^ = 0.028) 
and uniform stimulus ensemble, (b) Pinwheel density time courses for four different simulations (parameters 
as in a; gray traces, individual realizations; black trace, simulation in a; red trace, mean value) (c) OPMs 
and their power spectra in a representative simulation of Eq. ([3| with r(x) = and a/A = 0.12 (r^ = 0.57), 
other parameters as in a. (d) Pinwheel density time courses for four different simulations (parameters as in 
c; gray traces, individual realizations; black trace, simulation in c; red trace, mean value) 
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All together, the EN Eq. ([sj [4| and in particular the set of ground states of the EN model and their 
stability regions appear almost identical when considering either a circular or a uniform orientation 
stimulus ensemble. We found two different optima depending on the parameter regime, orientation stripes 
for larger interaction ranges and quadratic pinwheel crystals for shorter interaction ranges. In addition, the 
EN dynamics appears to be unchanged by the presence of unoriented stimuli. 



Appendix II 

Grid-like stimulus ensembles 

Refs. 64 65]) performed simulations with stimuli distributed in regular intervals in feature space, called 
grid-like ensemble. For comparison, we also performed deterministic annealing simulations with grid-like 
stimulus sets of varying size with non-periodic boundary conditions (see Methods). For these grid-like 
stimulus patterns, a competition between stripes and rhombs is observed. Notably, these are the only two 
stable states identified by our analysis for the circular stimulus ensemble. For non-periodic boundary 
conditions, rhombic pinwheel arrangements seem energetically favored for grid-like stimuli, almost 
independently of the size of the stimulus set. The average pinwheel density for N = 100 x 100 x 8 stimuli 
was p = 2.8. As expected from the predominantly rhombic arrangement of pinwheels, NN-pinwheel 
distances concentrate around half the typical column spacing and pinwheel pairs at short distances are not 



observed (Fig. 26 3). With these features, the maps obtained substantially differ from the experimentally 



observed pinwheel statistics |38] . 

The discrete EN model with 4th derivative 

In previous studies of the EN model, alternative definitions of the continuity term in the EN model have 
been explored |64j. A general continuity term for the spatially continuous formulation of the EN for 0PM 
and retinotopy is a linear differential operator which will suppress the emergence of high-frequencies during 
the EN dynamics. A finite-wavelength instability is expected in this case, although the precise expressions 
for the critical a and the typical wavelength will differ. As linear terms do not enter in the higher-order 
derivatives of the EN functional, changing the continuity term is not expected to alter the stability results 
obtained in the present study. 

To numerically test this expected robustness of our results for the EN model with discrete fixed sets of 



stimuli (see Figs. |18|and 19), we also performed simulations using deterministic annealing with a fourth 



derivative stencil (see Methods). Figure 27 illustrates that the results almost perfectly match the ones for 
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Figure 25: Phase diagram of the EN model for the joint mapping of visual space and orientation 
with a uniform orientation stimulus ensemble, (a) Regions of the ?7^-(j/A-plane in which n-ECPs or 
rPWCs have minimal energy, (b) Pinwheel density time courses for four different simulations of Eqs. ([sj |4| 
with r = 0.1, cr/A = 0.13(7^ = 0.51), r]r = V (grey traces, individual realizations; red trace, mean value; black 
trace, realization shown in c). (c) OPMs (upper row), their power spectra (middle row), and retinotopic 
maps (lower row) in a simulation of Eqs. (|3j|4|; parameters as in b. (d) Pinwheel density time courses for 
four different simulations of Eqs. ([sj [4| with r = 0.1, a/A = 0.3 (77 = 0.03), rjr = rj (grey traces, individual 
realizations; red trace, mean value; black trace, realization shown in e). (e) OPMs (upper row), their power 
spectra (middle row), and retinotopic maps (lower row) in a simulation of Eqs. ([3j|4|; parameters as in d. 
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Nonperiodic boundary conditions 

(grid-like ensemble) 

a N = 20x20x8 stinnuli p = 8 b 
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Figure 26: The EN model with deterministic anneahng and stimuli, distributed on a grid in 
feature space, (a) OPMs (left) and retinotopic maps (right) for A/" = 20x20x8 (upper row), A/" = 50x50x8 
(middle row) and N = 100 x 100 x 8 (lower row) stimuli and non-periodic boundary conditions (annealing 
rate x = 0.999). (3 is the continuity parameter in the conventional definition of the EN model (see Methods, 
Eq. (|46|) and is scaled, such that a comparable number of columns emerges in all simulation for each N. (b) 
Pinwheel densities of EN solutions for different numbers of stimuli, (annealing rate x = 0.999). Crosses mark 
individual simulations, red line indicates average values, (c) Statistics of nearest neighbor pinwheel distances 
for pinwheels of (upper left) arbitrary and (upper right) opposite and equal charge for 100 x 100 x 8 stimuli 
and non-periodic boundary conditions, averaged over four simulations (red curves). Black curves represent 
fits to the experimental data from 38 . Lower left: Standard deviations (SD) of pinwheel densities estimated 
from randomly selected regions in the 0PM. Black dashed curve indicates SD for a two-dimensional Poisson 
process of equal density. 
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the 2nd order derivative, considered in the main part of our paper (Figs. 18, 19 and 26 in the Appendix). 
When anneahng with periodic boundary conditions, the solutions very much resemble our gradient descent 
dynamics simulations with Laplacian term. The larger the set of stimuli, the more stripe-like are the 



orientation preference maps obtained (Fig. 27 i) and consequently pinwheel densities decrease (Fig. 27 3, 
upper right). The exponent for the SD is considerably lower than for the Poisson process (Fig. [27| 3, upper 
right). Typical NN-pinwheel distances concentrate around half the typical column spacing and in particular 



pinwheel pairs with short distances lack completely (Fig. 273, lower left and right). 

For non-periodic boundary conditions and random stimuli, we found that retinotopic distortions are much 
more pronounced. They however decreased with increasing number of stimuli. For large stimulus numbers, 
we observed stripe-like orientation preference domains which are interspersed with lattice-like pinwheel 
arrangements (see Fig. [27)3), lower row, upper left corner of the 0PM). Similarly to the periodic boundary 
conditions, short distance pinwheel pairs occur much less frequently than in the experimentally observed 
maps, indicating an increased regularity in the pinwheel arrangements compared to realistic OPMs (Fig. 



271, lower left and right). This regularity also manifests itself in a smaller exponent of the SD compared to 



the Poisson process (Fig. |27p). 



Simulations with grid-like stimulus as e.g. used in 64 65 displayed a strong tendency towards rhombic 



pinwheel arrangements analogous to the 2nd derivative case (Fig. [26^, f) 



Appendix III 

Strength of retinotopic coupling 

In our manuscript, we have shown that retinotopic distortions only have a weak influence on the optima of 



the EN model as well as its dynamics (see Figs. 10 and 12). Here, we quantify the influence of retinotopic 



distortions on the pattern formation process by comparing the angle-dependent interaction function for 
retinotopic coupling gr{(^) with angle-dependent interaction function of the EN model with flxed 
retinotopy. We use the ratio 

„ \\9r(-)h 



\\9i-)\\2 

as a measure to quantify the influence of retinotopic distortions. || • II2 denotes the 2-Norm, 

ll/(-)ll2= r f{a)da. 
Jo 

If c is larger than one, gr{o^) dominates the total interaction function gr{oi) + g{a) and retinotopic 
distortions may strongly influence the layout and stability of solutions of the EN model. On the other hand. 
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Figure 27: The EN model with deterministic anneahng and fourth derivative stencil, (a) OPMs 
(left) and retinotopic maps (right) for N = 10^ (upper row), N = 10^ (middle row) and N = 10^ (lower row) 
stimuli, non-periodic boundary conditions, and annealing rate x = 0.999. f3 is the continuity parameter in 
the conventional definition of the EN model (see Methods, Eq. (46)) and is scaled, such that a comparable 
number of columns is emerging in the simulations for each stimulus set. (b) Pinwheel densities (upper left) of 
EN solutions, standard deviations (SD) of pinwheel densities estimated from randomly selected regions in the 
solutions (upper right). Crosses mark individual simulations, red line indicates average values. Black dashed 
curve indicates SD for a two-dimensional Poisson process of equal density. Statistics of nearest neighbor 
pinwheel distances for pinwheels of arbitrary (lower left) and (lower right) opposite and equal charge for 
10^ random stimuli and periodic boundary conditions, averaged over four simulations (red curves). Black 
curves represent fits to the experimental data from [38j. (c) As a but for non-periodic boundary conditions, 
(d) As b, but for non-periodic boundary conditions, (e) OPMs (left) and retinotopic maps (right) for 
= 20 X 20 X 8 (upper row). A' = 50 x 50 x 8 (middle row) and A^ = 100 x 100 x 8 (lower row) stimuh, 
non-periodic boundary conditions, annealing rates x = 0.999. (f) As b, but for non-periodic boundary 
conditions and grid-like stimuli. 
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Figure 28: Strength of coupling between orientation map and retinotopy in the EN model, (a) 

Ratio of ||^r(-)ll2/lb(-)ll2 in the S4-cr/A-plane for the EN model at threshold and r] = rjr- (b) As a, but for 
r]r = 0, i.e. strongest coupling. Note the logarithmic scaling of the colormap. 



if c is small, the solutions and their stability properties are expected to not change substantially when 



including variable retinotopy into the EN model. Figure 28 displays the parameter c in the S4-cr/A-plane 



for the EN model at threshold for two different conditions, rj = rjr and rjr = 0. In the latter case, 
retinotopic distortions are expected to have the strongest impact. However, in both cases, c <C 1, in almost 
all of the parameter space, implying little influence of retinotopic deviations. Only for small <j/A and small 



54, c is larger than one. As shown in Fig. 12, this leads to slight deformations of the stability regions for 
rhombs, and stripes in this region of parameter space but does not result in novel optimal solutions. 
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